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Abstract 


An asymptotic technique is developed for analysing the 
propagation and dissipation of wave-like solutions to finite 
difference equations. It is shown that for each fixed 
complex frequency there are usually several wave solutions 
with different wavenumbers and the slowly varying amplitude 
of each satisfies an asymptotic amplitude equation which 
includes the effects of smoothly varying coefficients in the 
finite difference equation's. The local group velocity 
appears in this equation as the velocity of convection of 
the amplitude. Asymptotic boundary conditions coupling the 
amplitudes of the different wave solutions are also derived. 

A wavepacket theory is developed which predicts the 
motion, and interaction at boundaries, of wavepackets, 
wave-like disturbances of finite length. Comparison with 
numerical experiments demonstrates the success and 
limitations of the theory. 

Finally an asymptotic global stability analysis is 
developed which gives results which agree with other 
stability analyses and which can be applied to a wider range 
of problems . 
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1 . Introduction 


Consider the following very simple problem and 
numerical solution. The partial differential equation is 


3 u 
3 t 


+ 


c 


3 u 
3 x 


0 


where c is a positive constant. The domain 
considered is 0<x<1 . The initial condition is 


( 1 . 1 ) 


u(x,0) - exp [-200 ( x-0 . 5 ) 2 ] cos(kx) 


( 1 . 2 ) 


with k-80. This form of distribution is usually 
called a wavepacket. The cos(kx) term defines the 
oscillation of a group of waves and the exp [ - 2 0 0 ( x *0 . 5 ) 2 ] 
term is an amplitu'*'* 'envelope'. 


The upstream condition is 

u( 0 , t ) - 0 (1.3) 


The solution of this problem is 

u ( x-ct , 0 ) ct < x < 1 

0 0 < x < ct 


u(x , t ) 


(1.4) 


The numerical solution uses a uniform grid with 
computational domain 0<j<200 and a trapezoidal scheme 


,n+ 1 


U . - + T ( U . 


,n+ 1 


j + 1 


where 


+ 0 

c A t 
A x 


j + 1 


) - (U 




(1.5) 


( 1 . 6 ) 


In this example r»1. The initial condition is 


U 


u(x , 0) 


(1.7) 


and the upstream boundary condition is 


8 


U 


( 1 . 8 ) 


In addition a numerical boundary condition is 
required at the downstream boundary. For this condition 
space extrapolation is used. 


U 


200 


U 


199 


(1.9) 


Figure 1 shows the numerical solution at intervals 
of 60 time steps with each plot drawn to the same scale. 

The first two plots show the initial w?vepacket travelling 
downstream in the direction of the physical characteristic. 
Corresponding wavecrests are labelled a —e and it can be seen 
that the propagation velocity for the wave crests is greater 
than for the amplitude envelope. Note for example that the 
amplitude maximum lies approximately midway between crests b 
and d at n»60 but at n»120 the. maximum is clearly nearer 
crest b. At n*180 the numerical disturbance is interacting 
with the downstream boundary. The solution appears to be 
the sum of two waves, one with the original wavelength , and 
one with a very much shorter wavelength. At n*240 there is 
a reflected wavepacket of wavelength slightly greater than 
2. and the plots at n»300,360 show that this wavepacket 
travels back up the domain at approximately the same speed 
as the original wavepacket. This solution is clearly 
numerical and not physical since the analytic, physical 
solution moves from left to right across the domain and then 
out the downstream boundary. The analytic equation does not 
have any solutions with waves travelling from right to left. 
At n=420 the wavepacket is interacting with the upstream 
boundary, and at n-480 there is a reflected wavepacket with 
the original wavelength. This completes one cycle. If the 
solution was continued the wavepacket would travel down to 
the downstream boundary and then reflect again into a 
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n«60 Wavepocket Disturbance 
With Constant Frequency and 
Finite Amplitude Envelope. 

n«l20 The Amplitude Envelope 
Has Moved at a Slightly Smaller 
Velocity Than the Wavecrests. 

n « 180 The Wovepacket is 
Interacting With the Downstream 
Boundary. 


0 


n « 240 The Reflected Wovepacket 

Has Smaller Amplitude and a 

Wavelength Slightly Greater Than 2. 


0 




n * 300 The Wovepacket 
Travels Upstream. 


0 




n « 360 The Wovepacket 
Continues To Trovel Upstream. 




n* 420 The Wavepocket is 
Interacting With the Upstream 
Boundary. 

n * 480 The Reflected Wovepacket 
Has the Same Wavelength os 
the Original Wovepacket. 


0 


0 


.25 


.50 .75 1.00 

X 


FIGURE 1, NUMERICAL SOLUTION OF CONVECTION EQUATION 
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wavepacket with short wavelength, and a decreased amplitude, 
travelling upstream. 

The qualitative and quantitative prediction of the 
behaviour of numerical solutions in problems such as the 
above is one of the two objectives of this papar. The 
second objective is a global stability analysis 
incorporating boundary conditions and smoothly varying 
coef f wtients and predicting both stability and accurate 
asymptotic estimates of convergence rates. 

To achieve these aims a technique is developed to 
analyse the approximate time evolution of an amplitude 
modulated wave, i.e. a wave with fixed frequency and a 
slowly varying amplituae. Chapter 2 derives the theory for 
partial differential equations, while chapter 3 derives the 
theory for finite difference equations incorporating 
smoothly varying coefficients and boundary conditions. In 
the case of dispersive, non-dissipative wave propagation, it 
is found that the amplitude is convected at the local group 
velocity, a principle which is well understood in partial 
differential equations. 

Chapter 4 applies the theory to the motion of 
wavepackets which are wave-like disturbances of finite 
length and constant frequency such as in the earlier 
example. Chapters 5 and 5 derive global stability analyses 
with different levels of asymptotic approximation. Chapters 
7-9 develop further topics and examples including 
comparisons between numerical experiments and theoretical 
predictions . 

Throughout this paper a finite operator notation is 
used which greatly simplifies analysis and is a neccessity 
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for general proofs. Since there is no universally accepted 
standard notation Appendix 1 details the notation uaed. 

Very little previous work appears to have been done 
along the lines of thi3 paper. The concept of group velocity 
in partial differential equations is well understood and is 
explained in many texts The asymptotic approach of 

chapter 2 is not common due to the advantages of other 
methods but is discussed by Whitham 1 • Kentzer ^ has 
discussed the use of group velocity in analysing finite 
difference equations but does not derive a general equation 
for the amplitude or calculate the quantitative effects of 
boundary conditions. Vichnev e 1 3 ley and Bowles derive the 
the group velocity in finite difference equations using an 
approach which is valid only for constant coefficients. 

They also derive amplitude reflection coe f f f i ci e n t 3 at 
boundaries and discuss soma of the examples given in this 


paper. 
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2. Amplitude Analysis of Partial Differential Equations 


2.1 Fourier Analysis 


Consider a homogeneous partial differential equation 


L u ( x , t ) * 0 


-®<x<»,t>0 


( 2 . 1 ) 


where L is a constant linear differential operator 
defined by , 


L 


: f— 

' m n ' 3 x , 


m 


, - — 1 0 
l3tj 


( 2 . 2 ) 


ra , n 

and the 


C are constants, 
mn 


An e igenf unction of the operator I is a function 
u(x,t) satisf/ing 


L u ■ X u 


(2.3) 


where X is a constant called the eigen value . 

An ei genmod e is a solution cf the homogeneous 
equation (2.1) i.e. it is an eigenfunction with eigenvalue 
zero. 


exp[i(kx-ul. )] = ik exp [i( kx-ut) ] 


3 x 


(2.4a! 


— exp [ i ( kx-ut ) ] * -iu exp [ i( kx-ut ) ] 


(2.4b) 


3 m 3 n , m n 

— exp l i ( kx-u t ) ) - ( ik ) (-iu) exp [ i ( kx-ut ) ] (2.4c) 

i X i t 
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^HIG'NAL P«'Jc 15 
OF POOR QUALITY 


V fa ") m fa ) n 

) C — — exp i(kx-ut)] 

(_ mn 1.3 xj 1,3 tj 
m , n 


C (ik) m (-i(-!) n exp [i( kx-ut) ] 


mn 


( 2 . 4d) 


m , n 


Thus exp[i( kx-ut) ] is an eigenfunction of 

a a a^ m (a > l n 

- — , - — / t — 1 1 - — and L with eigenvalues ik , -iu , 

3x 3t [3xJ |3tJ 

(ik)®(-iu) n , and C (ik)®(-iu) n respectively 

L mn 

m , n 

Hence , 

u(x,t) * exp [ i ( kx-o t ) ] 

is an exact solution of (2.1) provided 

0 


( 2 . 5 ) 


y 


m , n 


C (ik)® (-iw) n 
mn 


( 2 . 6 ) 


This relation between k and 

u 

is 

called the 



dispersion relation. 






Examples of dispersion relations 

are ; 



Surface waves on deep water 

u 2 

- 

1 1 

( 2 . 

7a) 

Acoustic waves 

(j* 

= 

c 2 k 2 

( 2 . 

7b) 

Waves propagating along a waveguide 

u* 

9 

c 2 ( k 1 + kJ ) 

( 2 . 

7c) 

where g , c aid k 0 are constants. 






A general solution of (2.1) 

is 

a 

superposition 

of 


eigenmodes which in the case of a partial 

differential 



equation is expressed as an integral 

over 

all the 




wavenumbers k of the sum of all the eigen modes with 
wavenumber k. 


u ( x , t ) 


N 

A (k) exp [i(kx-u t)] dk 
/ n n 

n» 


( 2 . 8 ) 


I 
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If the dispersion relation is of order N in u , i.e. 
it contains powers of u up to u N , then 

(J x (k) , u 2 ( k ) , , u N ( k ) 

are the N values of u which satisfy the dispersion relation 
for a given value of k and the A^ik) are the corresponding 
constant amplitudes of those eigenmodes . 

If A^(k) is non-zero for all k,n then a neccessary and 
sufficient condition for u(x,t) to remain bounded and not 

increase exponentially is that each eigenmode must remain 
bounded. Splitting u into its real and imaginary components 
gives , 

u ■ u„ ♦ iu _ (2.9) 

R I 

exp [-iut] * exp[-iu t + u t] (2.10) 

R X 

Thus the condition that every eigenmode remain 
bounded, and hence a general solution remain bounded, is 

Uj < 0 for all k,n. 

This analysis is lacking in three respects. The 

first is that in some situations the initial d.sturbance is 

zero except for a finite region and one wants to know the 

time evolution of this disturbance, in particular the 

propagation velocity for the energy. The second failing is 

that when the initial-value problem is replaced by an 

initial-value / boundary-value problem with boundary 

conditions at x * 0 , 1 there is no easy way to include the 

effect of the boundary conditions in this stability analysis. 

The third failing is that exp ( i ( kx -u t ) ] is an eigenfunction 

of L only when the coefficients C are constant. The 

mn 

analysis breaks down when the coefficients are non-constant. 


r; 

15 


\ZJ 


The resolution of these problems requires the analysis of a 
wave of constant frequency with an amplitude which varies 
over a characteristic length scale much greater than the 
wav e length . 
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2.2 Asymptotic Amplitude Equation 


The problem now being considered is 

L(x) u(X/t) *0 -«<x<®,t>0 (2.11) 

where L(x) is a non-constant linear differential operator 
defined by. 


L ( x ) ■ ) C 

/ , mn 

m , n 

and the coefficients 
x . 


(X) 


(U 

m 

f 3 1 

WxJ 


Utj 


( 2 . 12 ) 


C (x) are slowly varying functions of 
mn 


The theory calculates the approximate evolution of a 
wavetrain with waves of a constant frequency u and a slowly 
varying amplitude, so u(x,t? is written as , 

u(x,t) » A(x,t) exp[i?(x,t)] (2.13) 


where A(x,t) is the slowly varying amplitude and 
? ( x , t ) is the phase of the wave which is related to the 
frequency u and wavenumber k by 
3T 

TT - -u (2.14) 


11 

3 x 


k 


The frequency u is constant but the wavenumber 
will vary slowly with x because of the slowly varying 


(2.15) 

k 


coefficients C (x) so the above relations can be integrated 
mn 

to giv‘e , 


/ 


* <x,t) 


k( K ) <35 - ut 


(2.16) 


jKIGINAL PACE i9 

OF POOR QUALITY 


To explain the asymptotic approximations which are 

made two characteristic length scales and and one 

characteristic time scale are defined. is the length 

scale for variations in k, is the length scale for 

variations in the amplitude A, and T a is the time scale for 

variations in A. Numerical values for L, , L. and T are 

k A A 

are given by , 


mm 


3 k 
3 x 


( 2 . 17a) 


L • min 
A 


3 A 
3 x 


(2.17b) 


mm 


3 A 
3 t 


(2.17c) 


The asymptotic approximations used in this theory are 


L, >> k~ l 
k 


L >> k" 1 
A 


T >> u~ l 
A 


which imply 


T 1 <K 

3 x 


3 A 


<< Ak 


3 A 


<< Au 


3 x 3 t 

A Taylor series expansion of A and ? about a point 


( x 0 , t „ ) gives 


3 A 

A ( x , t ) - A„ ♦ — (x-x„) + 

o X A 


3 A 
3 t, 


( t-t. ) + H . 0 . T 


(2.18) 


f (x,t) - f. - u( t-t. ) + k(5 ) d* 


/ 


- u ( t-t. ) + 


A 

/ 


3 k 


~ ( 5-x 0 ) + H.O.T ] at 

d x, 


[k 0 + 


18 


ORFG?^L F fS 
OF POOR UHY 

1 3 k 

- f 9 - u(t-t 0 ) + Jc 0 (x-x 0 ) + - — (x-x, l 1 + H.O.T. (2.19) 

2 3 X 0 

Subscript 9 denotes terms evaluated at (x 9 ,t 9 ). 

The H.O.T. , higher order terms, includes terras like 

\ — Y > TT7 » 7 — T A which are O { A ( L k ) “ 1 , A ( T u ) " * , A ( L, k ) ” * } 
3x3t3x A A k 

and are neglected in this asymptotic approximation. 


exp [if ( x , t ) ] ■ 

i 3 k 

exp|if„ + ik 9 (x-x 9 ) - iu(t-t 9 ) + - — (x-x 9 ) 2 + H.O.T 


* exp [ i f 9 + ik 9 (x-x 9 ) - iu(t-t 0 )] I 1 + -r (x-x 9 ) 

V, 2 o X 0 


. , f , i 3 k 
■-t a ) 1 1 + - 7 

+ H.O.T. 


* 1 


(2.20) 


Hence , 


u ( x , t ) » exp [ i f 9 +ik 9 ( x-x 9 ) -iu(t-t 9 )] 1 + — - — (x-x 9 ) 2 

2 o X ( 


i 3 k 


3 A , „ 3 A , 

a, * Tx 9 ( x-x ® * n 0 (t ' t * ) 


+ H.O.T. 


exp [if, + ik 9 (x-x 9 ) - iu(t-t 9 )J 


3 A 


3 A 


( t-t 9 ) + ^ A. (x-x-) 1 


3 x - (X “ X ° 5 + 3 t ' u_v - 9 ' ^ 2 3x, 


+ H.O.T. ( 2.21 ) 

To evaluate derivatives of u(x,t) at (x 9 ,t 9 ) a 
two-variable version of Leibnitz's rule is used. 


r 3 m f 3 

— ; ] — | ( f(x,t) g(x,t) ] 

m n 

Z l 

p*0 q»0 


m ! 


p ! ( m-p ) ! q ! ( n-q ) J 


f 


'aj p raj q 

,3xJ [3 tj 


<1 ^ m_P (3 > n-q 

3 x 


U t 
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; 

Of POOK QUALITY 


& 


Let 


f(x,t) ■ exp[if 0 + ik a ( x-x a ) - iu(t-t a )] 


9 ( x / 1 ) - A, + ~ (x-x 0 ) + (t-t,) + ^ A 0 ~~ (x-x,) 1 


at, 


2 0 3 x, 


Then , 


fa 'I 
l lax. 


atj 


( ik a ) P (-iu ) q exp [if, ] 


9 a * A 9 

a g ^ a a 
a x a ax, 


a_2 

at. 


a a 

at. 


111 

3x J 


. . ak 

1 • a x. 


(2.22) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 
(2.29) 


and all other derivatives of g(x,t) evaluated at (x 0 ,t 0 ) 
are zero. Hence , 


in 


— — | u ( x , t ) 

a x i a t 


exp [ if a ] • 


. ... , o , . n 3 A , . , .o—1 n 3 A ... ,o. . n — 1 

A „ ( lx # ) (-iu ) + m — ( i Jc 0 ) ( *iu ) + n (ik a ) (-iu) 

a X , 


m(m- 1 ) . 3 k , , . . m-2 

* 1 — — *• n. 


( -iu ) 


at, 


+ H. 0 . T . 


(2.30) 


and 


30 
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ORIGINAL F.' 

OF POOR QUALITY 


L ( x ) u( x , t ) - exp [ if ] 


e 

m , n 


mn 


m n 

A(ik) (-iu) + 


3 A ... .m-1, . .n . 3 A ... m. J .n-1 

m — ( ilc ) (-iu) + n — (ik) (*iu) 

. m( m- 1 ) . 3 k , . , , m- 2 , , .n ] .. _ _ 

i A — (ik) ( ”iu ) j + H.O.T. 

To satisfy the homogeneous equation (2.12) the 
amplitude A(x,t) must satisfy , 

3 A 3 A 

a 0 ( k , u , x ) A + a t (k,u,x) — + aj(k,u,x) — + a,(k,u,x) 


0 + rf . 0 . T . 


where , 
a # ( k , u , x ) 

a ! ( k,u , x ) 


a s ( k , u , x ) 


and a](k/U,x) 


E 

m , n 

E 

m , n 
3 a , 


C (x) ( ik ) “ ( - iu ) " 

mn 


C (x) n ( ik ) m ( -iu ) n 1 
mn 


3 u 


E 

m , n 


C (x) m (ik) m-1 (-iu )" 
mn 


- i 


3 a , 

Tk 


E 


# , . m( m- 1 ) , ,, .m-2, . , 

C (x) i (ik) (-iu) 

mn 2 . 


m , n 

i 3 1 a, 


2 3 k 2 

Because of the asymptotic assumptions , 


A » U-‘ || . 


k " 1 ^ , k * 1 A ^ 

3 x 3 x 


(2.31) 



(2.32) 
2 . 33a ) 

2.33b) 
2 . 33c) 
2.33d) 


so (2.32) can only be satisfied if 
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a, ( k,u ,x) - 0 


(2.34) 


This is the dispersion relation between k and u . 
k is now a slowly varying function of x due to the slow 
variation in the coefficients. Thus the characteristic 
length scale is related to some characteristic length 

scale L for variations in the coefficients. 

C 

Neglecting the H.O.T. and dividing by a ; gives the 
asymptotic amplitude equation. 


3 A 

at + c g 

3 A 

r— - e A 
3 x 

(2.35) 

where c ■ a. / 

g 

a i 

(2.36) 

and t • - a, 

3k . 

7x ' a ‘ 

(2.37) 


If c is real the left hand side of (2.35) is a 

g 

Lagrangian-type total time derivative with respect to an 


observer moving with velocity c 


If the coefficients C 


g mn 

am all constant, k is constant, e - 0 and so the amplitude 

A is constant along rays moving with velocity c^ . In a 

wavepacket individual wavecrests move with phase velocity 

u/k , usually denoted c , but the wavepacket, or amplitude 

P 

For this reason c is 

g 


'envelope' , moves with velocity c 


called the group velocity. 

Because the group velocity is the propagation 
velocity for the amplitude and energy of the wavepacket the 
group velocity is often more important than the phase 
velocity. One example is the Sommerfeld radiation condition 
which states that the waves generated by a fixed source have 
a group velocity directed away from the source. In some 
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unusual cases the phase velocity of the waves is actually 
directed towards the source. A second example is that the 
group velocity never exceeds the speed of light which is the 
limiting speed of propagation of information while the phase 
velocity can exceed the speed of light , as happens in wave 
propagation along an electromagnetic waveguide. 


To link this derivation of group velocity to other 
derivations the dispersion relation (2.34) is differentiated 
with x held constant. 


da 0 


3 a 9 
3 k 

0 


dk 



du 


Hence 


3 cj ) 3 a a 

k 3 kj x const 3k 



(2.30) 


* a i / a ; 


(2.39) 


The most common method of showing that 


UkJ 


x const 

is the group velocity uses the method of stationary phase 


which is well explained in the available literature [1,2] . 
The usual approach is to combine the dispersion relation , 
the .definition of the group velocity and some physical 
principle suc< as energy conservation to calculate the 
propagation of energy. The approach given above is not 
usually used partly because sometimes the exact partial 
differential equation is not known and the dispersion 
relation has been determined by asymptotic methods (e.g. 
water waves) or from empirical data (e.g. seismic waves). 
This approach is however suited to analysing finte 
difference schemes in which the exact finite difference 
equations are known and there is no general equivalent to 
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the principle of energy conservation. 


3. Amplitude Analysis of Finite Dif feranco Equations 


3.1 Fourier Analysis 

Consider a homogeneous finite difference equation 

L u" - 0 (3.1) 

As explained in the appendix A.1, L can always be 
expressed as a sum of step operators. 



m , p 

where the coefficients C are constants, but it is 

mp 

often more simply expressed as a polynomial of finite 
difference operators written symbol_cally as, 

L = (3.3) 

The eigenfunctions and eigenmodes of L are defined 
exactly as in 52.1. The finite operators all have the same 
eigenfunctions, exp [ i ( j • -nQ ) ] . * and Q are related to the 

wavenumber k and frequency u of the physical wave being 
modelled by, 


1 - kAx (3.4) 

0 ■ uit (3.5) 

As shown in the appendix A.1 

E^ exp [ i ( j * -nQ ) ] - exp(it) exp ( i ( j « -nfl ) ] (3.6a) 

exp [ i ( j • -nQ ) ] ■ 2i sin(#/2) exp [ i ( j * -nQ ) ] (3.6b) 

u exp [ i ( j s -nQ ) ] - cos($/2) exp [ i ( j $ -nQ ) ] (3.6c) 


E t exp [ i ( j ♦ -nQ ) ] - exp(-iQ) exp ( i ( j * -nQ ) ] 


( 3 .6d) 


V 
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5 t exp [ i ( j ♦ -nQ ) ] - -2i sin(Q/2) exp [ i ( j 9 -nQ ) ] (3.6e) 

u t exp [ i ( j 9 -nQ ) ] - cos ( Q /2 ) exp [ i ( j o -nQ ) ] (3.5f) 

E E exp [ i ( j * -nQ ) ] - exp [ i ( m* -pQ ) ] exp [ i ( j 9 -nQ ) ] (3. tig) 

nx p t 


Thus exp ( i ( j 9 -nQ ) ] is an eigenfunction of 


6 , u 
x x 


U and E E with eigenvalues 
t mx pt 


x x x t t 
exp(io), 2i sin(9/2), cos(i/2) , exp(-iQ), -2i sin(Q/2), 
cos ( 8/2 ) , exp [ i ( m$ -pQ ) ] respectively, and L has eigenvalue 


£ 

m , p 


C exp [ i ( m9 -pQ ) ] or 
m p 


P ( exp ( i 9 ) , 2 i sin ' 9 /2 ) , cos ( 9 /2 ) , 

exp(-iQ),-2i sin { Q/2 ) , cos ( Q/2 ) ] 
depending which expression for L is used. 


» exp [ i ( j 9 -nQ ) ] (3.7) 

is an exact solution of (3.1) provided 

C exp [ i ( m9 -pQ ) ] - 0 (3.8) 

m p 

m , p 

This is the dispersion relation between 9 and Q. 



Since exp[2niri] ■ 1 for all integers n, 9 + 2nu is 
equivalent to 9 so only solutions in the ranges 
- w < Re ( 9 ) < if 
-if < Re ( Q ) < t 

need be considered. 


If L involves P+ 1 time levelj and 
the dispersion relation is a polynomial of 
exp(-iQ) and of degree M in exp(i9). Thus 
tnere are P corresponding values of Q, and 


M+ 1 spatial nodes 
degree P in 
for & given 9 
for a given Q 


> 


I 
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there are M corresponding values of «. 

A general solution of (3.1J for periodic boudary 
conditions is a superposition of eigeniuodes. 

U n » ) 7 A <♦) exp[i(j*-nfl )] • 9 ) 

j L-t l_> P P 

• P" 1 

The $ summation is a summation over all the values 
of * which satisfy the periodic boundary conditions, and 
the p summation is over the P values of a corresponding to 
each value of * . 

If A (!) is non-zero for all «,p then a necessary and 
P 

sufficient condition fcr u” to remain bounded and not 

j 

increase exponentially is that each eigenmode must remain 
bounded. Splitting a into its real and imaginary components 
giv es , 

a - a ♦ in ( 3 . io ) 

R i 

exp[-inS] ■ exp[-inQ + na ] (3.11) 

K 1 

Thus the condition that every eigenmode remain 
bounded , and hence a general solution remain bounded , is 
Q < 0 for all » , p . 

This analysis is lacking in the same three respects 
as the analysis of partial differential equations by 
eigenmode expansion in the last chapter. The analysis 
gives no information about the movement of an initially 
localised disturbance, cannot incorporate boundary 
conditions or analyse schemes with non-constant 
coef ficients . 
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3.2 Asymptotic Amplitude Equation 


In this section the coefficients C in the 

mp 

definition of L (3.2) are assumed to be slowly varying 

functions of j. The analysis is performed in computational 
space with coordinates (j,n) in whicn the grid spacing is 
Aj-1, &n*1. Variations in mesh spacing in physical 
coordinates are incorporated directly into the variable 
coefficients of the finite difference equations. 


The theory calculates the approximate evolution of a 
wavetrain with waves of a constant frequency Q and a slowly 

varying amplitude , so o" is written as , 


U* - A(j,n) exp(if(j,n)] (3.12) 

where A(j,n) is the slowly varying amplitude and 
f(j,n) is the phase of the wave and is related to the 
frequency 3 and waven.aiber $ by 


11 

3 n 


3 


(3.13) 


11 
3 j 


t 


(3.14) 


which can be integrated to give, 

j 

?( j ,n) - *( 5 )dS-nQ (3.15) 

0 

As in §2.2 two characteristic length scales, L for 

$ 

variations in $ , and L for variations in the amplitude A 

A 

and one characteristic time scale T for variations in the 

A 

amplitude A , can be defined with numerical values beii.g 
given by. 


ORIGINAL PA\ 

OF POOR QLAlIV 


L ^ - min 


mm 


/ ii 
/ 3j 

1 / aj J 

t a - min ( A / f* J 

The asymptotic approximations are 


L 9 >> 1 


which imply 

»• « i 


l a >> 1 r A >> i 


3 A 

— << A 

3 D 


3 A 


<< A 


( 3 . 16a) 


(3.16b) 


(3.16c) 


3 j 3 j 3 n 

A Taylor series expansion of A ar.d ? about a point 


( j o ' n 0 ) gives, 


3 A 


3 A 


A(j 0 +m,n # +p) * A 0 + m — + p — + H.O.T 

d XI # 


3 j, 


(317) 


j 0 +m 


( j «, ♦m , n 9 *p ) * ? # - pfl + j 9 ( ? ) d£ 

3 o 


j« + o» 


f 0 - pQ + 


r 

j 


[♦# + |4 <5-j, ) + H.O.T ] d? 

d Jo 


- pQ + m« 8 +■ 2 il + H.O.T. (3.18) 

2 3 3 o 

Subscript 0 denotes terms evaluated at (j 0 #rt 0 ) . 

The H.O.T. , higher order terms, includes terms like 
3 2 A 3 2 A 3 2 » 


3j 2 3n 2 3j 2 


A which are 0 { AL a ” 2 , AT a " 2 , AL^ " 2 } and are 


neglected in this asymptotic approximation. 
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exp [ if ( j 0 +ra , n 0 +p ) ] 

im 2 3 ® 

■ exp [ if, - ipA + im® 0 + — — rr + H.O.T ] 

2 3 j , 


exp [ if „ - ip® + im® 0 ] 

Hence , 


1 + 


im 2 3 ® 
2 3 j, 


u n 0 + P » exp[if # ] exp [ i ( m® # -pA ) ] 
3 o +m 


1 + 


3 A 


3 A 


Ao + ® TT + P 7T 


3 j 


3 n. 


+ H.O.T. (3.19) 


im 2 3 ® 

2 3 j, 

+• H.O.T. 


* exp [ i 7 9 ] exp [ilni, -pA ) ] 


A, + m — - + p 

d J o 

+ H.O.T. 


3 A 3 A im 2 A 3_® 


3 n , 


Hence , 


L j U » exp [if] 


T. 

m , p 


C ( j ) exp [ i ( m® -pA ) ] • 

mp 


. 3 A 3 A in 2 . 3 « 

A + m — + p - — + — - A — 

3 j r 3n 2 3 ] 


+ H.O.T. 


3 3o 
(3.20) 


(3.21) 


To satisfy the homogeneous equation (3.1) the 
amplitude A(j,n) must satisfy. 


a 9 (®,n,j) A + a x (®,A,j) ~ + a,(®,A,j) + a,(®,A,j) A |~ 


3 3 

0 + H.O.T 


where , 

a 0 ( ® , A , j ) 

a x ( ® , A , j ) 


) C ( j ) exp [ i ( m® -pA ) ] 
L-i ™P 


m , p 


} C * ( j ) p exp [ i ( m® -pA ) ] 
!—> rap 


m , p 
i 


3_a o i 


3 A 


® , j const 


(3.22) 


( 3 .23a) 


( 3 .23b) 
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a 2 ( • , Q , j ) - 


and a , ( 9 , Q , j ) - 


In the above derivat:.on of a 9 ra^a^a, the general 
shift operator expression for (3.2) is used. In 

applications it is more convenient to use the finite 
operator polynomial expression (3.3). a, is obtained by 
replacing each operator with its corresponding eigenvalue 
and then a w a 2 ,a, are calculated by differentiating a 9 . 


) C (j) m exp [i(m*-pQ ) ] 
L . j ®P 


m , p 


-ifii.1 

* J 3 , j const 


( 3 .230 


C (j) — exp [ i ( mo -pa ) ] 
mp 2 


L 

m ,p 

2 la ♦ 1 ; a , j 


const 


(3.23d) 


Because of the asymptotic assumptions 


A >> 


a a 


a a 
3 j 


3 0 


a n a j a j 

so (3.22) can only be satisfied if 


a o < * * 2 » j ) 


+ 0 { L “ l ,L “ l ,T " l 
0 A A 


(3.24) 


This is the asymptotic form of the dispersion 
relation between o and a and will usually be satisfied by 

setting a 9 identically equal to zero. o is now a slwwly 


varying function of j due to the slow variation in the 

coefficients. The characteristic length scale L is related 

0 

to some characteristic length scale for variations in the 
coefficients. 

Neglecting the H.O.T. and dividing by a t gives the 
asymptotic amplitude equation. 
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3 A 

3 A 


3 n r 9 

3 j 

where 

r g * 

/ a 

and 

c - - ( 



e A 


ii 
3 j 


(3.25) 

(3.26) 

(3.27) 


Differentiating (3.24) with Q held constant gives, 
'3 a, 


da 


f 3 a 

•I <34 + — 0 dj 

' v 3 <p J Q,j const ^3j ) 4 , Q const 


- 0 ♦ H.O.T. 

Hence, neglecting the higher order terms, 

'3 a, 


(3.28) 


3 j 


f 3a ' 


3 j ) 4 ,Q const 


13# 


Q , j const 


i ( — *) 

^ 3 j J 4 , 8 const 


so 


-i a, 


f 3_a ( 

3 j J4 , 2 const 


a. a 


i “j 


(3.29) 

(3.30) 


If r is real the left hand side of (3.25) is a 
9 

Lagrangian-type total time derivative with respect to an 


observer moving with velocity r . Thus the amplitude A is 


being convected with velocity r^ in computational space 


Differentiating (3.24) with j held constant gives, 
(3 a, 


da , 


5 A ' 

, — -® I dB + -- 0 <34 

(.38 J 4,j const ^ 3 « J Q,j const 


- 0 + H.O.T. 

Hence, neglecting the higher order terms, 


(3.31) 


3_Q 

3 ip J j const 


3_a o 


3 4 


8 , j const 


(l*o) 

3 a 


4 , j const 


a, / a t 


■ r 


<£ 


(3.32) 
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Substituting for Q and p using (3.5a,b), 

9 (uAt) 
f g 3 ( k A x ) 


A__t 3 jj 
A x 3 k 


c At 

, 

A x 


(3.33) 


Thus r is the CFL number corresponding to the group 
9 

velocity in physical space of the propagating numerical 

wave. It is the number of spatial mesh intervals which a 
localised disturbance travels in one time step. For the 

rest of the paper r is called the group CFL number. 
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3 . 3 Examples 


The model problem which is considered is 

-« < x < ® 

Three different methods are analysed. 


3 u . . 3 u „ 

n + c(x> ■ 0 


(3.341 


3.3.1 Trapezoidal Scheme 


The trapezoidal scheme is 


L. 0 », * -li < 

At j j 4 A x 




I 1 U . . t-U . 

jl 5+1 3+ 

which can be written using operator notation as 


0 (3.35) 


1 

At 


j 

5 + <5 u 

t 2 2x t 


1 

3 


where the CFL number r is defined as 
c A t 


r . 
3 


Ax 


(3.36) 


(3.37) 


j 


and Ax. - - (x j+1 - x.^) 


(3.38) 


a„ is ootained by replacing the operators by their 
eigenvalues . 


so 


a, * -2i sin(Q/2) + ^ 2i sin($) cos ( Q/2 ) 

The dispersion relation is 
a a =0 

t a n ( 0 / 2 ) - ~ 9 i n ( • ) 

a. l , a, , a, are obtained by differentiating a 9 


( : . 39 ) 


(3.40) 

(3.41) 


34 


J? WGOK QwInUrt 



» cos(Q/2) + -j 3in(» ) sin(fl/2) 



* r cos ( <t ) cos ( 8/2 ) 



1 

2 3 « 

- ^ s in ( $ ) cos ( fl / 2 ) 
2 


Using the dispersion relation a x and a, can 
simplified . 

a l - cos(Q/2) + ^ sin($) sin(Q/2) 

- cos ( Q/2 ) + tan ( ft / 2 ) sin(fl/2) 

■ [ cos 1 (2/2 ) + sin J (Q/2) ] / cos(fl/2) 

- 1 / cos(Q/2) 


so 


and e 


a, * - - sin($) cos(fl/2) 

- - tan(Q/2) cos(fl/2) 
-* - sin ( r . /2 ) 

r g - 

* r cos ( $ ) cos 1 ( fl/2 ) 


-i a 


ii.) 


1 j v 3 j J $ , Q const 


/ a r a, ~ a fl / a t 


lr 


3 r 


sin(») cos(fl/2) [ — sin($) cos(S)/2) ]/[r 

d j 


(3.42) 


(3.43) 


(3.44) 

be 


(3.45) 


(3.46) 


(3.47) 


cos ( D ) ] 
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■ - 4 It sin*(f) cos*(8/2) / cos ( * ) (3.48) 

2 3 3 

There are three points of interest 
i ) sin ( tt ) - sin ( 9 ) 

so for all 8 there are two corresponding values of 
9 given by the dispersion relation, 

satisfying -n/2 < Re ( 9 x ) < »/2 
and *, -t - 9 x 

ii) For real 8 in the range 
0 < 8 < 2 tan" 1 ( r/2 ) 

0 < tan ( 8 /2 ) < r/2 
so 0 < sin ( ♦ ) < 1 

Thus and are both real and 

0 < 9 1 < t / 2 

so r ( ♦ . ) > 0 

g 1 

and ir/2 < 9 1 < n 
so r ( 9 2 ) < 0 

g 

Hence for every frequency in the given range there 
is one forward travelling wave, travelling in the same 
direction as the physical waves being modelled, and one 
backward travelling wave with wavelength less than 4Ax. 

iii) For real 8 in the range 
2 tan -1 ( r/2 ) < 8 < ir 


sin ( 9 ) > 1 

so and are complex 
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Let - »/2 + i* T 

Then sinf^ ) - cosh(« T ) 

so $ ^ is real and satisfies 

tan( a/2 ) - ^ cosh ( * J ) 

9 2 ■ w - ♦ 2 

■ ir / 2 - io ^ 

These are evanescent waves. If there are 
boundaries at j ■ 0,J and the boundary conditions force a 
steady oscillation with a frequency in the given range one 
wave will decay in amplitude exponentially away from the 
boundary at j m 0, while the other will decay exponentially 
away from the boundary at j-J. 


* 


j 



3.3.2 Box Scheme 
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OF POOR QUAUTt 


1 

At 


The box scheme is 


««r 


ii 

Ax 


* »r> 


- 0 


which may be written in operator notation as 

0 


It ‘ "x S t * r j“t J x 1 D jlt 


a«- -2i cos { 1/2 ) sin(Q/2) + 2ir cos(Q/2) sin(*/2) 
The dispersion relation is 


(3.49) 


(3.50) 

(3.51) 


tan(Q/2) ■ r tan(*/2) (3.52) 



i ^0 





a t ” 

an 





- 

cos ( * /2 ) 

cos ( Q / 2 ) + 

r sin ( n/2 ) 

sin ( 9 /2 ) 

(3.53) 


_ 1 3 a a 





a, ■ 
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- 

si n ( $ / 2 ) 

sin ( Q /2 ) + 

r cos ( Q/2 ) 

cos ( 4/2 ) 

(3.54) 

a, - 

1 3 2 a, 

2 aT 7 





- 

7 cos ( ♦ / 2 ) sin ( Q/2 ) 

4 

- 7 r cos ( Q/ 2 ) sin(*/2) 

4 


i 






' 8 *• 





- 

0 




(3.55) 

Thus r - a, 
9 1 

/ 4; 






. sin(»/2) s in ( Q / 2 ) + r cos(Q/2) coa(»/2) 
cos ( 9 / 2 ) co s ( Q / 2 ) ♦ r sin(Q/2) sin( */2 ) 
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tan ( 0/2 ) tan ( Q/2 ) » r 
1 ♦ r tan(*/2) tan(Q/2) 

1 tan* ( • /2 ) 

- r (3.56) 

1 + r* tan* ( ♦ /2 ) 

and c ■ 0 (3.57) 

There are two points of interest 

i) For each real value of Q there is one corresponding 

real value of • given by the dispersion relation and the 
group CFL number r is real and positive. 


ii) When r«1, r -r , so waves of all frequencies travel 

g 

at the same velocity as the physical waves being modelled. 
This is because when r» 1 the Box scheme reduces to. 


U 


n+ 1 

j + 1 



(3.58) 


which agrees exactly with the solution of the partial 
differential equation. 


u(x+ct,t) 


u ( X , 0 ) 


(3.59) 
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3.3.3 Backyard Eule: Scheme 


The backward Euler scheme is. 


— <u n+1 - u" ) * — 1 (o n+1 - o ntl ) - 0 

At j j 2 A x j+1 j-1 ' 0 

whxch may be written using operator notation as 


Sx > »r’ ■ ° 


a*" 1 - exp(ifl) + ir sin(a) 


The dispersion relation is 
exp(iQ] - 1 ■ ir sin(a) 

a - i li' 

1 3Q 

- exp(iQ) 


. 3 a, 

* * " 1 3a 


* r cos ( a ) 

1 . 3 2 a, 

a ’ " “ 2 1 77 ^ 


1 ii: 

2 3a' 


- s in ( a ) 


so r * a. 


- r cos ( a ) exp ( -ifl ) 


and . - -i >i3 conat / .... - a, / a. 


ir 3 r 

— sin(a) [ i — sin(a) ] / [ exp(ifl) r cos(a) 


1 i£ 

2 3 j 


3.60) 

3.61) 

3.62) 

3.63) 

3.64) 

3.65) 

3.66) 

3.67) 

1 


sin 2 (a) exp(-ifl) / cos ( a ) 


(3.68) 
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3.4 Asymptotic Boundary Conditions 


The general solution of 



is a sum of waves with different constant frequencies 9 and 
slowly varying wavenumber 9 and amplitude A 


0 


n 

j 



j 

exp [i( f v ( 5 ) d* -nfl ) ] 
J m 
0 


(3.69) 


The outer summation is over different values of 9, 
and the inner summation is over the M different values of $ 
which satisfy the dispersion relation for each 9. 


For each 9,m the amplitude A satisfies its 
asymptot _ amplitude equation on the interior of the 
computational eomain independent of all the other waves. 
All the waves of each frequency are however coupled by 
boundary conditions. 


Suppose a finite difference boundary condition at 


j * J is 


B u" - F n 

J 


(3.70) 


where B ;s a constant finite difference operator which can 
be expressed in operator polynomial form as 

B 5 WWW't' < 3 - 7 ” 


and F is a forcing function which can be expressed as a 
sum of inputs of different frequencies. 

r" - E ‘ ‘ 

9 


f ( 9 ) exp ( - i9 ) 


(3.72) 



OF POOR QUALITY 

Performing exactly the same asymptotic expansion as 
in the derivation of the asymptotic amplitude equation the 
boundary condition becomes, 

I t 

Q m» 1 
J 

exp [i( f $ ^ d£ -nfl ) ] * ) f ( fl ) exp(-inQ) + H.O.T. (3.73) 

0 Q 


. , 3 b 3 A . - _ . „ 

b A +1 — ; — m - i - — — : 


m 


3 Q 3 n 


3b 3 A i 

1 * — 7 — n* - — 

3 9 3 ] 2 


3 9 


3 2 b 


3 <p 2 m 3 j 


where b(fl,« ) * P T exp ( ±9 ),2i sin($ /2) ,cos(» /2), 

mum m m 


exp(-iJl),-2i sin ( fl/2 ), cos ( Cl/2 ) ] (3.74) 

The coefficients of exp(-inQ) in (3.73) must be 
asymptotically equal to zero for each Q so. 



3b 3 A 
i — — m 
3fl 3 n 


3b 3 A 
i 7 — — m 
3 4 3 j 


1 3 2 b , 3« I 

2 3 9 2 m3] 


J 

exo [i 9 d£ ] - f(Q) +H.O.T. (3.75) 

j m 

0 


This paper is primarily concerned with stability 
and convergence rates. When analysing perturbations from a 
steady state or constant amplitude oscillation the boundary 
condition f« r th perturbation has 

f ( Q ) - 0 (3.76) 


3ecause the zero order terms will usually dominate 
the normal form of the asymptotic boundary conditions is. 


A- J 

b ( 3 , 9 ) A ( J ) exp [ i f 


9 d? ] - 0 


(3.77) 
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The first order terms ~- m , , A q-^ 1 " are only 

d n d 3 d j 

important when, 

b(fl,«P ) - 0 ( T " l , L “ 1 , L ' 1 ) 

m A A $ 

As explained in §3.1 if the finite difference 
operator L in the interior scheme spans M+ 1 spatial levels 
there will be M values of $ given by the dispersion relation 
for a given value of fl. If the computational domain is 
0<j<J the interior scheme gives finite difference equations 
at J-M+1 nodes, so to complete the set of finite difference 
equations there must he M finite difference boundary 
conditions. Hence for each Q the asymptotic amplitude 
analysis gives M independent amplitude differential 
equations on the interior coupled at the boundaries by M 
boundary conditions. 
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3 . 5 Examples 


The same model problem as in 93.3 is considered, 


3 u 
3 t 


+ 



0 


(3.34) 


0 < x < x 

J 

c ( x ) > 0 

The analytic boundary condition is, 

u( 0 , t ) - F( t) (3.78) 

For perturbation analysis 

u( 0 ,t) - 0 (3.79) 

The finite difference scheme using the trapezoidal 
or backward Euler methods on the interior requires two 
finite difference boundary conditions. For perturbations 
the boundary condition at j = 0 is, 

ujj - 0 (3.80) 

The boundary condition at j*J is some form of 
extrapolation. Four of the most commonly used are analysed. 


3.5.1 Upstream Boundary 


U n - 0 (3.80) 

o 


B * 1 

so b * 1 


(3.81) 

(3.82) 


Hence Aj(0,n) + A,(0,n) » 0 


(3.83) 


In preparation for the theory developed in chapters 

4 and 5 it is useful to define R , the amplitude reflection 
coefficient as 
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< 3 / 


ORDINAL. 

or wqor q J/'*-* » ( 

(3.84) 


so in this example 

R Q - - 1 (3.85) 


A t ( 0 ,n) 
0 " Aj ( 0 , n ) 


3.5.2 Downstream Boundary : Space Extrapolation 


The space extrapolation boundary condition is 

u n ■ u n 

J J- 1 


(3.86) 


B - 1 - E 


-x 


E . 6 

- f x x 


(3.87) 


so b » 1 - exp ( -i$ ) 

= 2i exp(-i$/2) sin(*/2) 

Hence 


J 

r 


(3.88) 


siU* / 2 ) exp [ i , $ d£ - i * (J)/2 ] A (J,n) - 0 (3.89) 

m J m m m 

1 0 


The amplitude reflection coefficient R is defined 

U 


as , 


„ » * '•* ( j > n > 
rt J * Ai ( J,n) 

so in this example 
„ sin ( ♦ i / 2) 


sm ( 9 j / 2 ) 


exp 


(3.90) 


r j 

i J U; - * 2 ) d? - ■|(S l (J) - « 2 (J)) 

^ 0 

(3.91) 
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3.5.3 Downstream Boundary : Space-time Extrapolation 


£ 

r 

m* 1 


The space-time extrapolation boundary condition is 

(3.92) 


u n - u "" 1 

J J- 1 


B - 1 - E E 

-x -t 


b » 1 - exp ( - i$ + ifl) 

■ 2i exp [ i ( n ) / 2 ] sin [ ( | -fl ) / 2 ] 
Hence , 


(3.93) 


(3.94) 


sin 


* -Q', 


exp 


f J ' 

i ft d£ - i$ (J)/2 
j m hi 

k 0 


A ( J,n) - 0 (3.95) 

in 


sm 


l 2 


f J 


s m 


2 J 


exp 


j 


) dE. - -(•j ( J)-0, ( J) ) 


(3.96) 


3.5.4 Downstream Boundary ; Box Method 


In this example the Box method which was discussed 
in 53.3.2 as an interior scheme is now considered as a 
downstream boundary condition. 

- n 


1 “x 5 t + r “t S x 1 V* 


(3.97) 


B * E , E ■ [ u 5 + r u i 

-}-x -J-t x t t x 


(3.98) 


so b * exp [ i ( Q-* ) / 2 ] { -2i cos(t/2) sin(Q/2) + 

2 i r cos ( (1 / 2 ) sin(*/2) ) (3.99) 


Hence 
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L 


r sin ( t / 2 ) - tan(Q/2) cos($ / 2 ) 


m 


u 

exp[i f 9 dC - i ♦ ( J ) / 2 ] A ( J , n ) - 0 (3.100) 

J in tn m 


r sin( t i /2 i - tan(fl/2) cos($ 1 /2) . 

R J r sin ( * , / 2 ) - tan(Q/2) 003 ( 1 ) 2 / 2 ) 


exp 


y (#!-♦*) d? - ^(Oi 


( J)-« 2 ( J) ) 


(3.101) 
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4 . 

Ray Theory and Wavepacket-Particle 

Duality 


4_-_ 

1 Ray 

Theory 




In 

addition to the asymptotic 

approximations 

made 

in chapter 

3 

this chapter assumes that 

for all 

real 


wavenumbers $ 

, the frequency ft is real 

for all 

j and 

hence 

the group 

CFL 

number r is real. 

g 




r 

g 

-( 

an^ 
3 9J 

j const 



(3.32) 

A 

Lagrangian-type total time derivative in 


computational 

space is defined by. 




d 

dn 

= 

3 

3 n 

a 

+ r g 3j 



(4.1) 

so £2 

dn 

a 

ii 

3 n 

+ r |4 

g 3 3 





= 

r 

g 




(4.2) 

From 

the asymptotic amplitude 

equation 

(3.25) 

t 

dA 


3 A 

3 A 




dn 


3 n 

r 

g 3 3 





a 

e A 




(4.3) 

and using 

( 3 . 

26 ) and (3.29), 




d$ 


a ® 

a $ 




dn 


3 n 

g 3 3 







a ® 






g 

3 j 





- 

4 i 

iffiO - / a, 

(3 j J $ , ft const 2 





a 

i 1 

a L 1 

^3 j J * ,fl const 



(4.4) 
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A general initial value problem for a wave of 


frequency 4 

and 

wav enumber 

o (a, j) 

can 

be 

solved 

by 


integrating 

these equations 

(4.2)- 

( 4 . 

4 ) 

with 

initial 


conditions 










j ( 0 ) 

- 

j 0 







(4.5a) 

A ( 0 ) 

a 

A ( j , / 0 ) 







(4.5b) 

0(0) 

a 

0 ( Q # j , ) 







(4.5c) 

Each 

value of j generates 

a 

ray 

and 

all 

of 

the 


rays together cover the entire domain for n > 0. Figure 2 
shows the motion of some typical rays in computational 
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is a function of j so at a particular j all the 


rays have the same slope 


51 , 

dn 


Hence the time separation T of 


two rays, illustrated in figure 2 for rays 1 and 2, remains 


constant but the spatial separation L varies as r^ varies 


As explained in chapter 3 if the finite difference 
operator spans M+ 1 spatial nodes then for a particular 
value of Q there are M values of 9 which satisfy the 
dispersion relation. Define M+ to be the number of 
solutions 9 for which the group CFL number is positive, and 
similarly define M_ to be the number of solutions 9 for 
which the group CFL number is negative. Let the 
computational domain be 0 < j < J as usual. At j*J there 
are M + rays leaving the domain and M_ entering it. The 
amplitudes are related through the asymptotic boundary 
conditions each of which has the form. 



b( a, 9 ) 

m 


exp [ 



A ( J , n ) 

m 


0 


(3.37) 


Since the M + amplitudes of the rays leaving the 
domain are known and the M_ amplitudes of the rays 
entering the domain are unknown there must be boundary 
conditions to uniquely determine the amplitudes of the rays 
entering domain. Similarly at j*0 there mu3t be M + 
boundary conditions to uniquely determine the amplitudes of 
the rays entering the domain. 
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4.2 Wnvepacket- Particle Theory 


In terms of ray theory a wavepacket is a ray tube, a 
group of rays, along which the amplitude is non-zero. From 

the discussion in the last section the time length T. of the 

A 

wavepacket remains constant but its spatial length will 
vary whenever r has different values at the two ends of 

g 

the wavepacket. Provided << all the rays in the ray 
tube have approximately the same value for $(Q,j) so the 

motion of the wavepacket is given by. 


and 


dn 

d $ 
dn 


r ( Q , 9 ( j ) , j ) 

g 

i '3_a 9 'i 

a i [ 3 j j $ , Q const 


(4.2) 

(4.4) 


The energy, in physical space, of the wavepacket is 
defined as, 


E ( n ) - I | A ( x , t ) | 1 dx 

j n 


j I A( j , n ) | 2 


dx 

dj 


dj 


(4.6) 


Hence the wave energy density in computational 
space is detined to be, 


P ( j , n ) - a ( j ) | A < j , n ) | 2 


where a ( j ) - ~ 

dj 


(4.7) 

(4.8) 


Using (4.3) and the notation that A is the complex 
conjugate of A it follows that 

3 o 3 . 3 — 3 — 

— * tt ( r o) * t— ( aAA) + — (r a AA ) 

3n 3 j g 3n 3 j g 
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rs 
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— 3 A 3 A 3 A 3 A ~3 . 

- a A ( - — r — ) + a A ( r — + r ~) +• AA r a) 

3n g 3j 3n g 3j 3 j g 


- dA . dA - 3 . 

» aA — + aA — + AA — ( r a) 
dn dn 3 j g 


aAeA + aAeA + AA —t ( r a) 

3 j g 


t - 13 . 
e + e — ( r a ) 

a 3 D g 


Hence , 


dF f 3d 
dn J 3 n J 
0 


(4.9) 


- — ( r o ) + 

33 g 


' - 1 3 

e + c +■ — — ( r a ) 

a 33 g 


dj 


- r p 
g 


' 1 + ' * f r g° * I P dl 


(4.10) 


If the wav e packet is in the interior of the domain 
away from the boundaries the energy flux r^p at the 

boundaries j-0,J is zero. Also assuming as before that 

- 1 3 

L << L _ then e + e + — — (r •) ) is aporoximately constant 
AC a 3 ] g 

over the length of the wavepacket, so 

If ■ ( e * 7 * I h'v’) / ° dj 


1 3 , 

e + e ♦ - — ( r a ) 
a 3 j g 


(4.11) 


Thus equations (4.2), (4.., and (4.11) completely 

describe the motion of the wavepacket particle in the 
interior of the computational domain. 
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When the wavepacket reaches the boundary it 
interacts with the boundary conditions to produce one or 
more reflected wavepacket s with the same frequency but 
different wavenumbers. Tha only case for which it is easy 
to Incorporate boundary conditions is when M-2 and 

r ( *4 , j ) > 0 

r ( * , , j ) < 0 

g 

An example of a scheme satisfying this condition is 
the trapezoidal method which was applied to the model 
convection problem in the introduction. 


An additional assumption is that 
I r (fx #j) I # I r ( ♦, , j ) | << J 

so that it takes much more than one time step for a 
wavepacket to travel from one boundary to the other. 


Suppose that initially there is one wavepacket with 
wavenumber as in the introductory example. The 

wavepacket travels to the right with position and energy 


determined by the equations cf motion previously derived 
(4.2), (4.4) and (4.11). When the wavepacket reaches the 

boun !ary at j»*J a proportion of the energy is reflected 

into a left travelling wavepacket of frequency fl, wavenumber 
and energy Ej . Figure 1 illustrates this interaction. 


dE 

dn 


Equation (4.10) is 



2 

3 



3 ) 


\ 

J 


P d j 


(4.10) 


The outgoing energy flux is 
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r g ( • i # J ) o ; ( J , n ) '• r g (*!,J) a(j) | ( J , n ) | * (4.12) 

The incoming energy flux is 

r («.,J) o . ( J , n ) - r (*.,J) a(J) iA,(J,n)| l (4.13) 

g g 

The amplitude reflection coefficient R defined by 

J 


A 2 ( j,n) 

R ■ - - 

J A j ( J,n) 


(4.14) 


is a function of Q , $ l , <# 2 determined by the asymptotic 
boundary condition. 


The energy flux entering the reflected wavepacket S 

a factor 


r ( *, , J) 
9 

r ( 9 i , J ) 





greater, or less, than the energy flux leaving the incident 
wavepacket and so the total energy of the reflected 
wavepacket is given by. 


2 , 


r q < 9 , , J) 

r ( 9 j , 57 
g 



(4.15) 


The reflected wavepacket travels left according to 
the equations of motion for a wave 2 wavepacket until it 
reaches j-0 where i: is reflected into a right travelling 
wave 1 wavepacket. The reflected energy E t is given by, 


Ei * 


where 


r q ( tl ,0) 

I r q (* 2 , 0 ) 

A 2 V 0 , n ) 
R 0 " A, ( 0 ,n ) 



(4.16) 


(4.17) 


is determined by the asymptotic boundary condition. 



, 1 


t 

»• 




■) 


55 




ORIGINAL PAGE ,‘3 
0h p °OR QUALITY 


r ( • i , J ) 


( 9 l , J ) 


(4.15) 


& - r 


(♦, o) 


(4.2) 


d 0 I ^ _i f 3_a 0 'j 

dn 4 i k 3 j J « , fl 


const 


(4.4) 


e + e ■*■-—( r a ) Ej 
a 3 ] g J 2 


(4.11) 


r ( 9 . , 0 ) 

rr * 


r (♦, , 0 ) 
9 


R 0 I E ’ 


(4.16) 


In (a) and (c) the dispersion relation can be used 


The total number of time steps for a round trip from 


0 to J and back again to 0 is 


N = J [ r (»,,j) ] _1 - [ r ($,,j) ] " 1 dj 


(4.18) 


The energy growth of a wavepacket travelling from 0 


to J is given by 


d , „ . 1 dE, 

dn 4 E dn 


1 3 , 

e + e + - — ( r a ) 
a 33 g 


(4.19) 


d d di 

so -- In ( E ! ) - -- In ( E , ) / -r 1 
d] 4 dn * dn 


e + e 13, 

+ — ( r a) 

r r a 3 3 g 

9 9 
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e + e 


+ — [ In ( r a ) ] 

3 3 g 


(4.20) 


Hence , 


In [E l ( J) ] - In [ Ej ( 0 ) ] + y" dj + ln(r^a) 


so E x ( J) - 


r ( ♦ x , J) 

_2 i 

r (^ ,0) 


0 9 

J 


exp 


f *• 


( o ) 


(4.21) 


(4.22) 


Similarly the energy growth of the reflected 
wavepacket as it travels from J to 0 is 


so E 2 ( 0 ) = 


r ( * j , 0 ) 

_2 1 

r ( 0 2 , J ) 

g 


exp 


yj — 


E. ( J) (4.23) 


Combining ( 4 . 1 5 ) , ( 4 . 16 ) , ( 4 . 22 ) and (4.23) the round 
trip energy amplification factor X is 


r ( $ 2 , J ) 

_2 1 

r ( ♦ . , 0 ) 

g 


r ( <> . , 0 ) 

_2 !_ 

r^( $ 2 , J ) 


exp 

J 

/ 

( e + e ) 
r 


0 

l g J 

exp 

J 

f . 

f € + « 

J 

0 

1 r g 


dj 




r ( • . , J ) 

,_2 1 

r ( * , , J ) 

g 


r ( * . , 0 ) 

_i_J 

r ^ ( ♦ a , 0 ) 


0 J 


exp 


e + e 


g J 


e + e 


J. 


dj 


(4.24) 


where 


e + e 


l r 


g J 1 


is evaluated at , j 


and 


( e + e ^ 


is evaluated at <p 2 , j 


l g J* 

The condition for stability is 


X < 


1 
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The equivalent average decay rate a is 

a - - 1 In ( X ) (4.25) 

N 


When M> 2 and there are M + waves with positive r^ and 
M_ waves with negative r , one wavepacket with positive r 


reaching j*J produces M + reflected wavepackets with 
negative r , and one wavepacket with negative r read 
j=0 produces M+. reflected wavepackets with positive r 


Thus the total number of wavepackets increases with time 


exponentially. Since each wavepacket has finite lengtn 
this ultimately leads to the problem of determining the 
effect of interference between overlapping wavepackets. In 
general the sum of the energies of two wavepackets is not 
equal to the energy of the sum of the wavepackets. If the 
wavepackets are identical the latter is twice the former 
while if the second wavepacket has opposite sign to the 
first the latter is zero. 


58 


5. Asymptotic Stability and Convergence Analysis I 

5 . 1 Theory 

In addition to the asymptotic approximations in 
chapter 3 this chapter assumes that M* 2 an^ there is one 
boundary condition at both j=0 and j=J and if is real at 

j=0 then » ! and $ 2 are real over the whole domain. 

Examples of methods satisfying these conditions are 
the trapezoidal method applied to the model convective 
problem with variable CFL number r, and the backward Euler 
method applied to the model convective problem with constant 
CFL number. Methods which do not satisfy these conditions 
include Lax-Wendroff and Runge-Kutta type schemes. For 
these methods the general stability analysis of chapter 6 is 
require d . 

As explained in ?3.1 a standard Fourier series 
analysis of the problem with constant coefficients and 
periodic boundarv conditions shows that the eigenf requencies 
are 3(#) wnere $ is a real wavenumber satisfying the 
periodic boundary conditions and 3($) is the corresponding 
frequency giv en by the dispersion relation. 

The common use of a Fourier series analysis to 
predict the stability of problems with non-periodic boundary 
conditions implicitly assumes that for real wavenumber $, 

3 ( $ ) is a close approximation to an actual eigenf requency . 
This chapter follows that assumption, calculates a 
correction 3' to this 3($) due to the boundary conditions, 
and then determines the validity of the assumption for this 
particular class of methods based on the asymptotic errors. 
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For real $ 2 , $ 2 and the corresponding complex ft the 
general solution is , 



f 

n 1 

] 


f 

ft 1 

\ 

* A 1 ( j , n ) exp 

i 

\ d? -n 3 


+ A s ( j , n ) exp 

i 

/ 0, d£ -nB 




, 0 J 

J 



lo J 

J 


(5.1) 


where the amplitudes A t and A 2 satisfy the asymptotic 
amplitude equations 


3A- ( , 3A- 

7— 31 + ( r ) 7 — m * e A 

3n g m 3 j mm 


and boundary conditions 


f f 

b; ( a , * x ) exp 1 i J * 2 dj 


m- 1 , 2 


+ b , ( Q , 0 j ) exp 


* b i 

(3 ,« 2 ) A, ( 


1 


dj 

A t ( J , n ) 


J 



r 

J ^ 

r 

exp 

i 

1 •, dj 


L 

0 j 


(5.2) 


(5.3) 


Aj ( J,n) 


(5.4) 


A linear system of finite difference equations with 
time-independent coefficients and boundary conditions has 
eigenmode solutions of the form. 


U - . * ( U. ) . exp ( -inB ) 

.3 9 j 3 

where 3 is the comp 1 ex eigenf r equency and U is 

3 3 

time-independent . 

Suppose the frequency 3 in (5.1) is close to 3 
Define 3 1 by 


(5.5) 


3 „ = 3 + 3 

8 


(5.6) 


Thus equating (5.1) and (5.5) 


3 


. > UAU * i 


(U.). exp(-infl') - OF r'J^' ' 

o j 

' j 1 [ j 

A ,. ( j # n ) exp i I d? + A, ( j , n ) exp j i J $ , d£ (5.7) 

l 0 J < 0 

so A|(j,n) * exp(-inQ') A t ( j , 0 ) (5.8) 

A j ( j , n ) = exp(-inQ') A, ( j , 0 ) (5.9) 


(Ug). » A x ( j , 0 ) exp j i J $ l d£ + A, ( j , 0 ) exp i J ♦, d£ 

l 0 J l 0 


Substituting (5. 8), (5. 9) into (5.2) gives 


3 f t + ifl ' 'I 

— (A m (j,0)] - A (j,o) 

3 1 m r m 

l g Jm 

which can be integrated to give 


m = 1 , 2 


A ( J , 0 ) » A (0,0) exp 
si si 


e + iQ 


0 v 9 ) m J 


The boundary conditions (5. 3), (5. 4) then become in 


matrix form 


B A =» 0 


where A * 


| A, (0,0) 

Ui ( 0 ,0 ) 


r 3 1 1 a i 2 

i b j i B n 


B u - b. ( Q , $ t ) 


(5.15b) 
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(5.15c) 


I J 

i j dj 


f J 


B,, ■ b 2 ( fl , $ , ) exp 


♦ j dj 


exp 


exp 



e + iQ 


0 
( J 


dj 


K o K g 


( 5 . 1 5d ) 


( 5 . 1 5 e ) 


A non-zero solution A of (5.13) exists if, and only 


if. 


exp 


de t 3 

Hence 

J 


(5.16) 


iQ 


' / t r g (9 x ,3) ]’ 1 - [ r g U,,j) J" 1 dj 


b; ( ^ 2 ^ ( 11 f $ t ) 

b 2 ( 2 # ♦ i ) b r ( Q , $ 2 ) 


exp i 


i \ «, - 9 x dj 


r J 

r re y re y 

e * P | i I r ] - ,7 i 

{ o ^ <3 j i ^ qji 


(5.17) 

the right hand 3ide of (5.17) can be expressed as 
its magnitude 


bj(Q,9 ? ) b v ( Q > 1> t ) 
b,(Q,9 ; ) b j ( Q , 9 2 ) 


exp Re 


(i)-r-n 

' a/ l 


V o v qJ l ^ q j i j 

multiplied by a phase factor exp(i'f) where 

J f J N 

b, ( a , t, ) b , ( Q , 9 , ) f r,'c> r e ^ 

) b < a ., ) * i « ♦ *■ / ? * b « 

0 l o MM 

(5.18) 

If 9 X , <2 are chosen so that ?/2ir is an integer then 


(5.17) 


reduces 


to 
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fl ' 


i 

N 


In 


b, (flj, ) b, ( fl , * , ) 
b 2 ( Q , ♦ i ) b 2 ( Q , $ 2 ) 



(5.19) 


where as defined in chapter 4, 

J 

N - f [ r ]" 1 - [ r U 2 ,j) ]“‘ dj (4.18) 

J 9 9 

0 


The stability criterion is 

Im( Q a ) - Im( fl + Q 1 ) <0 (5.20) 

8 


Thus the frequency fl resulting from a normal Von 
Neumann analysis is corrected by an amount O' due to 
boundary conditions and variable coefficients. This 

approach, using 0 as an initial approximation to fl the 

9 

actual eigenf requency , is valid provided the asymptotic 
errors are small compared to O'. 


The asymptotic error is 0( L^ -2 , T^“ 2 ) * 0( L c ” 2 , ft' 2 ) 


Now N - 0 — 

r 

g 

so if r < < J then N >> 1 

g 

and hence fl' << 1 except near frequencies for which 

b 2 \fl,? 2 ) bj ( fl > $ i ) 
b, ( Q , * 2 ) b 2 ( Q , $ j , 

is zero, or infinite, which usually occurs at fl«0 . However 
these frequencies are heavily damped by the boundary 
conditions and so an accurate estimate of their 
eigenf requencies is not essential. This method gives 
accurate asymptotic values near the critical frequencies 
which are least damped and which therefore determine the 


m 


3/ 
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overall spectral radius of 


the scheme . 



5 . 2 Example 
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This example is the Backward Euler method applied to 
the model convective problem with constant CFL number r and 
space extrapolation at the downstream boundary. 

The dispersion relation is 


exp(ifl) - 1 ■ ir sin($) 
so if is real, $ 2 ■ ir - $ 2 is also real 


(3.63) 


b l (fl, » ) 


(3.82) 


bj(Q,$) * 2i exp(-i«/2) sin(*/2) 


(3.88) 


b j ( 0 i t ] ) b t (0,$; ) 

b, ( a , • x ) b t ( fl , ♦ , ) 


exp [-i ( o , t ) / 2 ] 


sin ( o , /2 ) 
sin ( $ i /2 ) 


- exp [-i( ) / 2 ] 


cos((, / 2 ) 
sin ( o x /2 ) 


- axp [ -i ( i , -« x ) /2 ] cot(« l /2) (5.21) 


c -0 

- ♦ l )/2 + J(l, - «i) 

1/2) (ir - 2% x ) (5.22) 

f / 2 tt - n, where n is an integer, implies 
• i - t/ 2 ~ 2irn / ( 2J- 1 ) (5.23) 

The group CFL numbers are 
r^(«j) ■ r cos( »! ) exp(-ifl) 

* r cos(«, ) / ( 1 + ir sin( «, ) ) 


Since — - 0 , 

3 3 

Hence f - ( $ 2 

* (J - 
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Q ' 


so 


Im( Q 


Thus 


so 


» r cos)*! ) ( 1 - ir sinO, ) ) / (1 + 

■ - r cosl^ )( 1 - ir sin(0i))/(1 + 
Hence , 

ir cost*, ) (1 - ir sin(*, )) 

■ - A . - . ■! ]_ n 

2 J ( 1 + r l sin J (* v ) ) 

exp[-2 Im(fl)] - | exp(i3) | 2 

* | 1 + ir sin ( * l / 2 ) 

* 1 + r J sin 1 ( o x /2 ) 


r 1 sin* ( o l ) ) 
r 1 s in* ( o t ) ) 

! cot ( o i /2 )j 


Im(Q) - - j In [ 1 + r*sin , (o l /2) ] 


(5.24) 

(5.25) 


(5.26) 


(5.27) 

(5.28) 


Hence 


O' ) - - - In [ 1+ r I sin J (« l /2) 


r cost », /2 ) 


In 


2 J ( 1 + r 1 sin 1 ( o x /2 ) 
Define the decay rate j to be 


- In( Q + O' ) 


For small o << 1 
a * r 1 o 1 + In ( 2/0 ) 


do 

do 

da 
d o 


r # * 2J * 


- i 


* 0 at o * ( 2r J ) "i" 


min p(o) * + f- In ( 8 r J ) 

4 J 4 J 

The spectral radius X is 


cot(0^/2) ] (5.29) 

(5.30) 

(5.31) 

(5.32) 

(5.33) 

(5.34) 


X 


max exp ( -<J ) 
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* max ( 1 -a ) 

» 1 - ~ [ 1 ♦ In ( 8 r J ) ] (5.35) 

4 J 

56.2 continues this example proving rigorously that 
every eigenmode is stable, and deriving an asymptotic 

expression for the decay rate of the eigenmodes near ^-0, 
showing as expected that the decay rate is greater than the 
minimum decay rate obtained above and so the above analysis 
is valid in calculating the spectral radius. 



' 13 

Or'pOOB QUAUfV 


6. Asymptotic Stability and Convergence Analysis II 


6.1 The or' 


This chapter continues the analysis of chapter 3 
without any additional assumptions or approximations. The 
eigenmodes of a linear system of finite difference equations 
with time- in depen dent coefficients and boundary conditions 
vary exponentially with t.'ime »o a general eigenmode can be 
written as 


M r j > 

0 n - exp(-inQ) ) A (j) exp i f * d£ 

, m J m 

m» 1 l 0 > 


( 6 . 1 ) 


Note that the amplitudes A are independent of n. 

si 

The time evolution of the eigenraode is contained solely in 

the term exp(-inQ). The complex amplitudes A each satisfy 

m 

their asymptotic amplitude equation. 


r - c A 

g 3 j m 


( 6 . 2 ) 


' J 

V J » ’ V 01 ** p / (r ] 


(6.3) 


j — 1 denotes — evaluated at * 
r r m J 

' oJ m q 


Now that A (J) is related to A (0) the M asymptotic 
m m 

boundary condions can be written in tensor summation form as 


B, A * 0 
km m 


(6.4) 


where A "A ( 0 ) 
m m 


(6.5) 


km 


b (Q,* ) if the k th boundary condition is applied 

k m 

at j - 0 

( 6 . 6 ) 


b.(Q,* ) «xp 


f J 

/I 

r e ' 

r 

V. 0 

^ 'J' 


m 


+ it dj 
m 


if the k th boundary condition is applied 

at j - J 


r 

The term b, (3,« ) exp / it 

k m J m 

V 0 


dj 


comes from (3.77) 



r J 

\ 


r 

ft y 

and the term exp 1 

J 

7 dj 


l 0 

t gj m j 


is the factor relating A (J) 

m 


to A ( 0 ) in (6.3). 

m 


A non-zero solution to (6.4) exists if, and only if, 
det ( B ) - 0 (6.7) 


Since all the elements of B are implicit functions 
of 3 this is the equation which determines the 
eigenfrequencies . 

If the coefficients are constant all of the A are 

m 

constant and so there are no asymptotic errors. The theory 

is then exact and is identical to the P-stability analysis 
of Beam, Warming and Yee (5]. If the coefficients are 

variable the asymptotic error is of order 0(L _l ). 

C C 

The 0{L "*) comes from neglecting second derivatives of A 
C m 

and # in the asymptotic amplitude equations. The 0(J" l L„" 1 ) 

comes from neglecting the first derivatives of A and * in 

m m 

the asymptotic boundary conditions. 


For all but the very simplest problems it will be 
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impossible to solve (6.7) analytically to obtain the 
eigenf requencies . Three possible numerical approaches are 
outlined below. 

6.1.1 Iterative Solution 


Suppose a n is a good approximation to an 
eigenf requency 3 and 
_n+ 1 


_ n . _ n 

a + Aa 


( 6 . 8 ) 


is to be a better approximation. 


The terms in the definition of B which change most 

Km 


rapidly with variations in a are exp 


i$ dj since these 
m 


are oscillatory functions because $ usually has a real 

m 


component of order 0(1). 

3« 

3Q [3djm 

■ (r ) " l 
\ g/ m 

n+ 1 n r \ - 1 

Hence 9 a 9 + (r A3 

m m \ g/ m 


(6.9) 


( 6 . 10 ) 


Thus 


Subscriot m means evaluated at $ ,j 

m 

Superscript n means evaluated at 3 n 


f J t 

f n + 1 

exo I / i$ d] 

i j a i 

a 

i 

ex p 

r 1 

f l *« dj 

f J 

expi j iA3 (i 

1 

g-> dj 

' 0 J 



l 0 J 

lo 

J 


( 6 . 11 ) 


Define 
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b, ( fl , $ ) if the k boundary condition 

km 

is applied at j-0 


m 

r j 

1 


f 

J 


. . _ n n . 

b. ( Q ,5 ) exp 

k m 

M 

lo 

i ^ a n _ , 

+ i5 d j 
m 

exp 

i A n 

/ ('.] 

in “ 1 
l m 


m J 



o ^ j 

■ / 


if the k th boundary condition is at j*J (6.12) 


Then choose AO to be the smallest root of 

det ( B n ) -0 (6.13) 

The method in chapter 5 performs one seep in this 
iterative procedure. As explained in chapter 5 the 

asymptotic error remaining after a correction AO is 
0 ( AQ" 1 ^^" 1 ). For constant coefficients this gives 

quadratic convergence to the true eigenf requency. For 

variable coefficients it will converge to a frequency which 

differs from the true eigenfrequency by an asymptotic error 

of order 0 ( L * 1 ) . 

C 

The solution procedure will fail as described in 

chapter 5 near frequencies for which 

the fractional variation in b, (Q,$ ) 

k m 

may exceed, the fractional variation 


b, ( 0 , $ ) is zero since 
k m 

is comparable to, or 
f J 


i n exp 


/* 

l 0 


15 dj 
m 


.1 


6.1.? Newton- Raph son Solution 


The Hewton-Raphson solution procedure is 

g n + 1 ^ ^ n de t [ B ( Q ° ) ) 

fr det [B( Q n ) ] 


(6.14) 


Although this seens straightforward there will be 
problems in practice because det(B)*Q has many roots and so 
det[B(fl n )] has many zeros near which the Newton- Ra ph son 
iterative procedure is badly behaved. 


6.1.3 Stability Domnin Method 
If z is defined as 

z = exp ( iQ ) (6.15) 

then the condition for stability is that det(B) has no zeros 
in | z S < 1. In very simple cases this criterion can be 

tested analytically (an example is given in 36.2). In more 
complicated cases because det[B(z)] is an analytic function 
of z the ’Principle of the Argument' method outlined in 
appendix A. 2 can be used to find whether det[B( z ) ] has any 
zeros in |z| < 1. To find the spectral radius X the same 
method can be used to find the largest X for which there are 
no zeros in |z| < X. 


6.2 Example 0 p POOR QUAUTT 

The example is the Trapezoidal method applied to 
the model convective problem with constant CFL number r 
space extrapolation at the downstream boundary. 

The dispersion relation is 



tan ( Q/ 2 ) 

r 

m — 

2 

s in ( $ ) 

( 3 

Define 

< » exp ( i$ ) 


( 6 


z * exp ( - 

in ) 


( 6 

Then 

tan ( 0/2 ) 

i 

exp(in/ 2 ) - exp(-iQ/ 2 ) 


i 

exp(in/ 2 ) + exp(-in/ 2 ) 




2 

1 - exp ( -in ) 




i 

1 + exp ( - in ) 




- 

2 1 -z 
i 1 + z 

(6 


s in ( $ ) » 

1 

2 i 

[ exp ( i* ) - exp ( -i« ) ] 



- 

1 

2 i 

( < ~ 1 ) 

( 6 


Thus the 

dispersion relation becomes 



1 - z r 

— =* - ( < -K 

1 + Z 4 

) 

( 6 


Note that 

if 

is one solution then tc 2 » -< 

- 1 

l 

the other solution. 

The condition for stability is 

that 

| z | < 

1 for all eigenmodes . 


Let z 

« R exp ( i 9 ) 


R , 9 real, R> 0 

( 6 

Then 

1-z 1 - 

1+z * 1 + 

R exp ( i 9 ) 
R exp ( i 9 ) 



= [1 - 

R 

exp( 18 )] [1 + R exp(-i9)] 



C 1 + 

R 

exp(i9)] [1 + R exp(-i9)] 



and 

.41 ) 

. 16 ) 
. 17 ) 

. 18 ) 

• 19 ) 

. 20 ) 
is 

.21 ) 
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1 - 2 i R sin ( 9 ) - R » 

1 + 2 R cos(3) + R* 

(1 - R 2 ) - 2i R sin ( 9 ) 

[1 + R cos ( 9 )] 1 + [R sin(9)J* 


( 6 . 22 ) 


Hence R < 1 


Re 


1 - 2 
1 + z 


> 0 


R > 1 => Re 4^ < 0 

1 + z 

Thus if r is positive then R 5 | z| < 1 if, and only 

<■■> Pj ] > 0 


if. Re 4— > 0 

1 + z 


Now let k ■ R exp(i9) (6.23) 

< - k" 1 * R cos ( 9 ) + iR sin(9) - R“ l cos(9) - iR sin(9) 


(R - R“ * ) cos ( 9 ) + i ( R + R" 1 ) s in ( 9 ) 


(6.24) 


so 


I 2 I < 1 


Re [ k “ 1 ] > 0 


either R > 1 

or R < 1 


-ir/2 < 9 < ir / 2 
ir/2 < 9 < 3ir / 2 


Figure 4 shows the stability domains for the 
different variables. 


The upstream boundary condition has 

b l * 1 (3.82) 


The downstream space extrapolation has 


Hence B = 


■ 1 - exp( -ifl ) (3.88) 

1 1 ' 
exp ( i J$ ) [1-exp(-io l )] exp ( i J$ ) [1-exp(-i« l )] 


1 


1 


, J (1 


( -<" 1 ) J ( 1 + <) J 


<~ l ) 


(6.25) 


V 
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det ( B ) = 0 implies 

k 2 J " ^ ( k - 1) - (-i) J (< + 1) » 0 

k+1 , J 2 J- 1 

7 - ( - 1 ) ic 


(6.26) 

(6.27) 


Consider the twc i nsw* .-ility regions I and II marked 
in figure 4 iii ) . 


so 


In region I 
i > 

< - 1 


K +> 1 


, . J 2 J- 1 

( - 1 ) < 


< - 1 


and | < | > 1 


so 


In region II | < + 1 | < | < 


< + 1 
K ~ 1 


, , . J 2J-1 
( - 1 ) « 


1 | and | k | < 1 


Hence there are no possible solutions of (6.32) in 
the two regions of instability. Every eigenmode is stable 
and so the method is stable. This example is a special 
example taken from a more general result, proved by Beam, 


Warming and Yee [5] , that the class of A-stable Beam-Warming 
multistep methods with q t ^ order space extrapolation at the 


downstream boundary is stable. 


Continuing this example asymptotic decay rates for 
small can be derived. 


$ ! a 0 * > < a 1 

Let < * 1 + <5 (6.29) 

Equation (6.31) becomes 

(1 + j) 2,1 " 1 ^ - (-1) J (2 ♦ 6 ) * 0 (6.29) 

Now ( 1 + ^) J a exp(x) (6.30) 

j-*® j 


so provided J >> 1 >> 5 
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Hence 


so 


Since 
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( 1 ♦ 5 ) 2 J - exp ( 2 Jd ) (6.31) 

exp ( 2 Jd ) * 2 ( - 1 ) J d ” 1 (2 + <5 ) ( 1 + 5) 

» 2 (-1 ) J d _l (6.32 ) 

If J is even this has a solution for real 5 
2 J 6 ■ ln( 2/6 ) (6.33) 


h lai2/,) 



±_ In 

4 J 



2 J 

In ( 2/d ) 



In ( J ) 
2 J 

+ terms of order 0( 

In [ In ( J ) ] 
J 

(6.34) 


- exp ( i* ! ) 

(6.35) 

* -id 


i In ( J ) 

3 • ' — 

2 J 


(6.36) 


Linearising the dispersion relation about #»u gives 
a r 9 ! 


ir In ( J ) 

3 _ -■■■■’ 

2 J 

The asymptotic decay rate is 

r In ( J ) 

a * 

2 J 

If J is odd the smallest d solutions are 
. r In ( J ) _ ir i 

o * ■ n — 

2 J 2 J 


(6.37) 


(6.38) 


(6.39) 


and the decay 


rate 


is the same 


7. Assorted Examples and Further Developments 


7.1 Instability of Backward Euler with Spacetime Extrapolation 

Consider the Backward Euler method applied to the 
model convective problem with constant CFL number r on 
domain 0 < j < J with space-time extrapolation at the 
downstream boundary and J >> r >> 1. 

The dispersion relation is 

exp(iQ) - 1 » ir sin(« ) (3.63) 

The upstream boundary has 

b ! ■ 1 (3.82) 

The downstream boundary with space-time 
extrapolation has 

bj = 1 - exp[i(Q-«)J (3.94) 

The eigenf r equency equation, det(B)-0 reduces to 

exp(iJ^ x ){1 - ex*. [ i ( Q -« t ) ] } * exp(iJ$ 2 ){1 - exp [ i ( 3 -« 2 ) ) } 

(7.1) 

Because t, * this can be written as 

exp(2iJ* l ) (1 - exp [ i ( 0 -« 2 ) ] } - (-1) J (1 + exp [ i ( a + ® 2 ) ] } * 0 

(7.2) 

This and the dispersion relation form two equations 
in the two unknowns Q and . Considering only the case in 
which J is eve one eigenmode is given by, 

& * ir + fl ' | Q ' | < » (7.3) 

so ir sin(9 l ) ■ -2 - ifl 1 0 ( S2 1 1 ) (7.4) 
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— - - ' + H.O.T. 

r r 


(7.5) 


Substituting into (7.3) with 8'*0 gives 
exp(-4J/r) (1 + exp(2/r)} - 1 + exp(-2/r) » 0 + H.O.T. (7.6) 

Expanding the exponentials under the approximation 
r >> 1 this reduces to. 


J--ln(r) + H.O.T. 
4 


(7.7) 


For a particular even value of J the value of r 
satisfying (7.7) makes this eigenmode asymptotically 
neutrally stable. To find whether increasing r makes it 
unstable or not (7.2) is differentiated by r. 
d D , 

2iJ — 1 exp(2iJi> l ) (1 - expli(J) - ^ )] ) 

+ exp ( 2 i J® , ) (i -j— x - i 7 ^} exp(i(fl - ®.)] 

ar ar 


- (i 4^ 1 + 1 exp [ i ( Q + ® i ) ] 

dr dr 1 


(7.8) 


This is evaluated at fl'-O, J * — ln(r) so 

4 


, , , d®. 1 , d®. dfl , , d® . dfl , 

dr r dr dr dr dr 


0 +H.O.T. (7.9) 


d fl d® . 

Hence — * ln(r) — 1 + H.O.T. 

dr dr 


2 i In ( r ) 


H.O.T 


(7.10) 


As r increases from the neutrally stable value, 
Im(fl) becomes negative so the eigenmode becomes unstable. 
Thus this eigenmode is stable only if 

7 In ( r ) < J 

4 

In numerical experiments it is found that this is 


© 


t 


the least stable eigenmode and so this condition is both 
necessary and sufficient. The condition is asymptotically 
equal to an exact stability condition derived by Beam, 
Warming and Yee [5] . 
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7.2 Optimum CFL Number For Trapezoidal Method 


Consider the Trapezoidal method applied to the 
model convective problem with constant CFL number r on 
domain 0 < j < J with space extrapolation at the downstream 
boundary. 

The dispersion relation is 

tan(Q/2) ■ ~ sin($) (3.52) 

The upstream boundary condition has 

b x » 1 (3.82) 


The downstream boundary has 

b, - 2i axp(-i$/2) sin($/2) (3.38) 

Following the stability analysis of chapter 5, e-0 

since r is constant, and the group CFL numbers are given by, 

r - r cos(<#.) cosMtf/2) (3.47) 

9 i 

r ■ r cos(«, ) cos 1 ( 3/2) 

<? i 

■ - r cos(* x ) cos 1 (3/2) (7.11) 


s o 


arg b. (Q ,», ) b x (Q ,$, ) 


J 

♦ dj 

0 



- ( J-1/2 ) ( •, - ) 


- (J-1/2) ( ir - 2*! ) 


(5.18) 


(7.12) 


?/2ir ■ n wheis n is an integer implies 
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a ' - 


where 


t 


i 

N 


i 

N 

i 

N 


N 


•» 

2* n 

2 

2 J- 1 

r 


In 

b, ( Q , ♦ , ) b , ( fl , * , ) 

b, ( Q , « , ) bj ( n , • , ) 




In [ sin($,/2) / sin(# l /2) ] 


In [ cot ( ♦ ! /2 ) ] 

J 

- i [ ( * l , j ) ] “ l - [r g ( * 
0 

2J 

r cos ( * t ) cos 1 ( Q/2 ) 


(7.13) 



f J 


♦ Re 

f (? ) * (1 

l 0 V <? ; i V 



J) 


(7.14) 

,j)]" 1 dj (5.20) 

(7.15) 


The decay rate a was defined as 
a - - im( a ♦ a • ) 


o ■ cos( ) cos 1 ( 3 / 2 ) In [cotl^/2) ] 

2 J 

r cosUi ) In [ cotU, / 2) ] 

2J ( 1 ♦ tan 1 ( 0/2 ) ] 

2 r cosd, ) In [ cot(» t /2) ] 

J [ 4 ♦ r 2 sin 1 ( # t ) ] 


(5.30) 


(7.16) 


Figure 5 shows the dacay rate j as a function of r 
for various values of $ . 


When a 3ystem of finite difference equations is 
solved using a time-independent approach to a steady state 
solution, usually the initial conditions and the final 
solution are smooth so the initial error is smooth i.e. the 
error is predominantly in the low wavenumber, long 
wavelength eigenmodes. For this problem suppose that the 


\ 






3/ 
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initial error is in the wavenumber range 0 < f < $ # where $ 0 
is a constant, $ 0 < if/2 


3n 

U 



sin(l ) 

4 + r J sin J ( « ) 


2r : sin(» ) cos 2 ( t ) 
[4 + r 2 sin 2 ( t ) ] 1 „ 


ln[cot($/2 ) ] 


cos ( $ ) 

sin ( $ ; ] 4 ♦ r 2 sin 2 (♦ )] 


(7.17) 


so over the range given above a has a minimum at <p , . To 
maximise the overall rate of convergence the CFL number r 


can be changed by altering At. 

3 o _ cos ( 3 ) [ 4 - r 2 sin 2 ( <t ) ] 
3r 2 J [ 4 + r 2 si n 2 ( « ) ] 

> 0 


In ( cot ( # /2 ) ] 


for 0 < r < 2/ sin ( $ ) 

* for 2/sin (♦ ) < r 

so a ( <p 0 ) is maximised by choosing 

2 

r * — : — 

sin(», ) 


(7.18) 


(7.19) 


The smoother the initial error the lower the value 
of 4 „ and the higher the optimum CFL number. If the initial 
error is not at all smooth with t 9 approaching ir/2, the 
optimum CFL number drops down towards 2. 
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7.3 Discontinuous CFL Number In Trapezoidal Method 


Consider the following Trapezoidal method applied to 
the model convective problem. 


[ «. 


+ lizk 

2 


r . 


U 7 + 

t x 2 


Ilk 


u n +* - 


u A ] u" 1 - 0 

t x j 


(7.20) 


Suppose that r is discontinuous 
r , j < 0 


r . 
3 


j > 0 


On the two sides of the discon u-.ouity a general 
solution of frequency fl can be written as, 

r A t ( j , n ) exp [ i ( j ® l -nQ ) ] + A* ( j , n ) exp [ i ( j $ , -nQ ) ] j < 0 

u“ * ^ (7.21) 

i Aj( j,n)exp[i( j»,-nQ)l + ( j , n ) exp [ i ( j 0 , -nQ ) ] j > 0 


®. and « 2 satisfy the dispersion relation 
r 

tan(Q/2) » — “ sin(®) (3.52) 

with ® j = ir - ® : 
and 0 < Re ( ® L ) < ir/2 

®j ar. d ®„ satisfy the dispersion relation 
r 

t an ( Q / 2 ) = -+ sin(e) (3.52) 

with 

and 0<Re(®j)<ir/2 
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The amplitudes Aj , A, are related to A, , A % at j»0 
by two equations . The first comes from the requirement that 
the two expressions for 0^ are equivalent at j*0. 


Thus Aj + A* * A, + A, 

The second requirement is that 

[6 + 4 r . , u 7 +4 r. . W A 1 U*? + 2 

1 t 2 t x 2 3 +J. t x 0 


(7.22) 


(7.23) 


which implies 

[5 + r u 7 ] 

t t x 0 


n-fl 

-[5 + r u A ] U_ 2 

t + t x 0 


(7.24) 


The left hand side involves U at j*0,-1 so the 

first expression for U n is used. The right hand side 
n 3 

involves at j»0,1 so the second expression is used. 
Neglecting derivatives of the amplitudes the resulting 

equation is , 

A v { -2i sin(Q/2) + r _cos ( 8 / 2 ) [ 1 - exp ( - i $ l ) ] > 

+ A, { -2i sin(fl/2) + r_cos ( 0/ 2 ) [1 -exp ( - i • , ) ] > 

* - A, { -2i sin(Q/2) + r ^cos ( n / 2 ) [ 1 - exp ( - i o , ) ] } (7.25) 

- A* { -2i sin(Q/2) + r + co s ( (2 / 2 ) [ 1 - exp ( - i $ , ) ] } 

Substituting for tan(£l/2) using the dispersion 
relations and using 


exp(i«) = cos($) + i sin( » ) 


(7.26) 


this reduces to 

r_{ Aj [1-co3( ) ] + A, [1-cos(», )] } 

» r + { A, [ 1-cos( I j )] + A* [ 1-cos ( ♦ , ) ] } 


(7.27 ) 
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In the general problem the two equations relate the 
amplitudes on either side of the discontinuity at j*0. Now 
consider the particular problem in which # a , , are 

all real and A, is zero. This is the situation when a 
wavepacket with wavenumber $ x , and positive group CFL number 
travels from j<0 to the interface at j-0 producing a 

transmitted wavepacket with wavenumber and positive group 

CFL number and a reflected wavepacket with wavenumber $ x and 

negative group CFL number. The reflection coefficient R and 

transmission coefficient T are defined by, 

F - < t / A x (7.28) 

T ■ A, / A x (7.29) 


Since A„»0 


A i + A 2 


(7.30) 


so r_ { A x [ 1-cos (#,) ] + A, [ 1-cos ( ♦ j ) ] } * r + [ A x +A 2 ] [ 1 -cos ( $ , ) ] 

(7.31 ) 

r_ [ 1 -cos ( $ x ) ] - r^tl-cosU, )] 


Hence A 


ind 


r_[ 1-cosUj ) ] - r f [1-cos(« s )] 

1 * R 

r_ [ 1-cos ( « 2 ) 1 - r _ l 1 -cos ( 9 x ) ] 

r_ [ 1-cos ( ft 2 ) 1 " r_ [ 1 -cos ( ft , ) ] 


(7.32) 


(7.33) 
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8 . 1 Degeneracy 


So far in this paper it has been implicitly assumed 
that no two eicenmodes have the same frequency 8 and 
wavenumber 8. This section considers the degenerate case 
in which this happens. An example is the Trapezoidal 
method. The dispersion relation is 


( 3.41 ) 
These are 


tan (8/2) - — sinU) 

which has wavenumber solutions $ and $ * ir-9 

iii 

ident. cal when 

9 X ' = - v/2 (3.1) 

which occurs at the degenerate frequency 8' given by, 

tan( 8/2 ) « t (8.2) 


r 

2 


. r a 8 v 

Note that r * j - — I 

g ^ 3 9 j j const 

* r cosU') cos 2 (Q'/2) 

* 0 (8.3) 

This is characteristic of degeneracy because in the 

neighbourhood of the double zero 


8-8 


so 


3 8 
3 $ 

2 a ( 
0 


a ( 9 - 9 ' ) 


- 9 ' ) 
at 9 * 9 ' 


constant 


(8.4) 


(8.5) 


For 8 8 ' the general solution for constant r is 
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■ A l exp [i ( j* ! -nfl ) ] + A, exp [ i ( j o , -nfl ) ] 

Consider the limit as fl approaches O' 


0-0 


^ fn - n ' 

' a ± I 


*■ 


SO 0 ! “O' » - ( 0 j - 0 ' ) 


( 8 . 6 ) 

(8.7) 

( 8 . 8 ) 


Hence A l exp[ij<> l ] + A t [ i j $ 2 1 

- exp [ i j $ 1 ] { Ajexptijt^- 0')] + A, exp [ i j ( ♦ , - $ ' ) ] > 
» expfijo'] { A l exp[ij(« I - $ ' ) ] + A,exp[ij($'- x ) ] > 
exp[ijo'] { ( A t +A 8 ) cos [ j ( o t -o ' ) ] + i ( A x -A, ) sir j ( o x -o ' ) ] } 


» exp [i j 0 ' ] ( A, + A* j } 


(8.9) 


where A, » A, + A. 


( 8 . 10 ) 


and A* = i(A x - A 2 ) ( 0 x - O') 


( 8 . 11 ) 


If A, and A* are now considered to be constants in 
the limiting process the general form of the degenerate 

eigenfunction is, 
n 


U" * (A, + A, j ) exp [ i ( jo ' -nQ ' ) ] 


( 8 . 12 ) 


Another way of deriving this result is through the 
asymptotic amplitude equation. With the amplitude A 

time-independent, r constant, and r equal to zero the 

9 

largest term is the second order derivative of A so, 


3 : A 

3 j 2 


(8.13) 


which implies A = A, + A* j (8.14) 

In all the examples I have analysed I have not yet 


found a degenerate eigenfunction which satisfies the 
asymptotic boundary conditions and so is a degenerate 
eigenmode. The degenerate frequency does however satisfy 
the determinant condition 

det(B) - 0 (6.7) 

because two of the wavenumbers are equal so their columns 
are identical and B is singular. This is not a problem in 
the stability analysis because the separation of 
eigenf requencie s , using the result from chapter 5 , is 

approximately, 


2ir 

N 


2ir 


J 

[ [r ( t x , j ) ] * l 


> 

[r ( * j ) 1 " l dj 


^ 0 J 

* 0 as Q * Q ' 

so there is usually a true eigenf requency which differs 
from the degenerate frequency by less than the overall 
asymptotic error. 


1 
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8.2 Eigenmodes and Eigenf requencie3 


In chapters 5 and 6 it was stated that the 
eigenmodes of a linear system of finite difference 
equations with time-independent coefficients and boundary 


conditions vary exponentially 
outlines the proof. 

Let U n - z n V. 

D j 

If the domain is 0 < j 
finite difference equations so 

C V « 0 


with time. This section 

(8.15) 

< J then there must be J+ 1 
V ^ satisfy 

(8.16) 


where V is the J+ 1 vector of elements and C is a ( J+ 1 ) * 

matrix whose elements are polynomials of z. 


For there to be a non-zero solution requires 

det ( 0-0 (8.17) 


This is the equation that determines the eigenvalues 
z of the eigenmodes. Apart from the problem of possible 
degeneracy the on j.y remaining difficulty is to show that the 
number of eigenmodes equals the number of independent 
initial conditions needed to start a numerical solution. 

For the Backward Euler method with space extrapolation at 
the downstream boundary there are J-1 independent initial 
conditions since 


J- 1 


( 8 . IS ) 


and 


(8.19) 


II u 
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r , z/ 2 

z- 1 r , z/2 


" r j-1 2/2 2 -' r j-i 2/2 

1 - 1 

/ 

( 8 . 20 ) 


so det(C) is a polynomial in z of degree J- 1 giving J- 1 

eigenmodes. Thus any solution can be expressed as a sum of 
eigennodes which vary exponentially with time. 


0 

fl 


► 


8.3 Other Asymptotic Approaches 


One approach which can be used when the variation 
in coefficients i a extremely small is to let, 

u" ■ A ( j , n ) exp [ i ( j« -nft ) ] (8.21) 

where ft and $ are both constant and satisfy 

a, (ft,* , j) - 0 (8.22) 

at some point j . 

The asymptotic amplitude equation is then 

. 3 A 3 A . _ . _ . 

a 9 A + a l - — *• a 2 — * 0 (8.23) 

d n d j 

where a 0 , a L and a, are defined and calculated as before. 


L 


A 


min 


I A 


/ 


3 A 'I 

3 3 J 


* 0 (a, / a 0 ) 


(8.24) 


so the fractional error using this method is O ( a # /a, ) 2 . 

If the variations in the coefficients are small this is fine 

but if the variations are 0(1) the fractional error is 0(1) 
i.e. the method fails to give accurate asymptotic 
approxima tioi.s . The asymptotic approach used in this paper 
allows total variations in the coefficients of 0(1) and only 
requires that the length scale of those variations is much 
greater than 1 . 


Another approach is to set 

u!J ■ A(j,n) exp t if ( j , n ) ] (8.25) 


with 


11 

S3 


■ * 


(8.26) 


and A, £2, 9 are all real and slowly varying. This approach 

is used extensively in the analysis of water waves and other 

partial differential equations with dispersion and very 

little dissipation. This approach applied to finite 

difference equations would give poor results because 

dissipation over one time step is 0(1) so if £2 is real A 

reflects this dissipation and so T * 0(1). The method 

A 

presented in this paper is able to use constant complex £2 

rather than variable real £2 as in the above method because 
the eigenmodes have constant complex ei gen f r equen c i e s 
provided the finite difference equations are 
time-independent . 


S. Wavepacket Teat Program 


9.1 Program Description 

The program solves the model convective problem 
3 u 3 u 


3 t 


c - — - 0 

d X 


(9.1) 


using a choice of Box or Trapezoidal methods on the domain, 

0 < j < 200 (9.2) 

and time step range, 

0 < n < 400 (9.3) 


The CFL number 


j 


c ,4t 

Ax . 
3 


(9.4) 


is specified by the user at j-0,200 and the program 
interpolates for intermediate values by fitting an 

exponential curve through r Q and r . 


r( j ) 


r 0 ex P 


' ^ 

^ 2 Oo" 1Mr 2OO /r 0 ) 


(9.5) 


Methods 1-3 are different Trapezoidal methods which 
are identical if r is constant. For these methods the 
program offers a choirs of space extrapolation, space-time 
extrapolation or box condition as the downstream boundary 
condition. Method 4 is the Box method. The wavepacket 
theory for eacli of these methods is derived in the next 
section . 


The upstream boundary condition is, 

ujj - 0 (9.6) 


The initial conditions are given by, 


2 / 
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U 


j 


A ( j , n ) exp [ i? ( j , n ) ] 


where A(j,n) - ox p [ - ( j - 1 0 0 ) * / 2 0 0 ] 
and t(j,n) - ) * i ( Q # j ) 


(9.7) 

(9.8) 

(9.9) 


This produces a wavepacket about 40 mesh points 
long. The user specifies 0 by specifying tne critical CFL 
number which is defined as, 

2 t a n ( Q / 2 ) (9.10) 


so 


when. 


cr 1 1 
fl - 2 tan* 1 (r 


/ 2 ) 


(9.11) 


cr i t 

The significance of the critical CFL number is that 


cri t 

the group CFL number for the Trapezoidal methods is zero 
Wave propagation can only occur for r > r crit; • 


(9.12) 


At each time step the program calculates 
'experimental' values for the position X(n) and energy S(n) 
of the wavepacket. If U was continuous X and E would be 
de fined aw, 


E( n) 


200 

f 

J 


I U(x,t ) | 2 dx 
n 


200 

f 

j 

0 


ff « 


I U( j ,n ) | 1 


(9.13) 
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'200 


and X ( n ) 


— • r i 

E(n) J 1 


U ( x , t ) | 2 x dx 
n 


0 

200 


iuT / i0(3 ' n)|1 X|J) If 


dj 


(9.14) 


so since U is discrete. X(n) and E(n) are defined as, 


E( 


2Q_0 

- 1 

j-0 


(9.15) 


X ( n ) - 


200 

— r i u n 

E( n) /. 1 j 


x , a 


j j 


where ^ ^ , > 

In the program c is taken to be constant so 
variations in r are due to variations in Ax. 


Ax 


c A t 


(9.16) 


(9.17) 


(9.18) 


so 


X , . - X , . 

3+1 3-1 




A physical domain 0<x<1 is used so, 

3 / 200 

, 5 , ('>•»)'■ 


(9.19) 


After completing the 400 time steps the program 
calculates predicted values for X(n) and E(n) using the 
wavepacket equations derived in the next section with the 
experimental values at n- 1 as initial cor ^itxons . 
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9.2 Theory 


Method 1 is a Trapezoidal method which Is second 
order accurate in Ax only when Ax is constant, 
r . 


[6 + — u 6 ] U n + ^ 

t 2 t 2x j 


r . 
3 


2 c A t 


X . -X . 

3-1 


(9.20) 


(9.21) 


Following the analyses of chapters 3 and 4 
a, * - 2 i s in ( 2 / 2 ) + ir cos(2/2) sin(t ) 


(9.22) 


a. * co3(!)/2) + sin(2/2) sinU) 


r cos ( Q/2 ) cos ( t 


(9.23) 

(9.24) 


- — cos ( fl/2 ) sin ( $ ) 
2 


(9.25) 


The dispersion relation is 


t a n ( 2 / 2 ) » - si.n(» ) 


so a 


(9.26) 

(9.27) 


and a! * cos(2/2) + sin(2/2) tan(2/2> 

- ( cos 2 ( Q/2 ) + sin* ( 0/ 2 ) ) / cos(2/2) 


* 

sec ( 2/ 2 ) 

(9.28) 

r ■ 

g 

a, / a t 

(3.26) 

a 

r cos ( $ ) cos 2 (2/2) 

(9.29) 

d i# 

dn 

i f — •! / . 

L 3 j J 2 , <# const 1 

(4.4) 


( -P- 


98 


ORIGINAL PAGE Rj 
OF POOR QUALITY 


V3/ 


e 



sin ( ♦ ) 


cos* (Q/ 2 ) 



2 > 9 


const 


/ 


a 




a, / a t 


■ - — 7 ^ cos 1 ( 0 / 2 ) 
2 3: 


sin ( * ) tan ( ♦ ) 


Since a is 


proportional to r~ l , 


1 3a 1 3 r 

a 3 j r 9j 


so 


1 — 
a 3 j 


( r a ) 
9 


r 


3 j 


(r /r) 

g 


(9 
( 3 
(9 


( 9 


* - r sin!; ! cos 1 ( 0 / 2 ) 
» - r sin($) cos J (Q/ 2 ) 

* cos 2 ( 0 / 2 ) sin(«) 


3 t 

3 j 



tan ( $ ) 


/ »: 

(9 


Hence the equations of motion for the wavepacket 
real 0 and 9 are , 


*1 

dn 

d$ 

dn 

d 

dn 


* r cos ( $ ) cos 2 ( Q/ 2 ) 

(9 

- - sin( * ) cos 2 (Q/ 2 ) 

3 3 

( 9 

ln(E) - 0 

(9 


Method 2 is also a Trapezoidal method which is 
second order accurate in Ax only when Ax is constant. 


3-f . 

[ 5 t * -J— H t 4 x * — 


n + £ 
0 . 2 
J 


30 ) 

30 ) 

31 ) 

32 ) 


33 ) 
for 

. 34 ) 
.30 ) 
• 35 ) 


2 


2 


0 


(9.36) 
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cA t 


X . -X . 

: + i 3 


1 f 

- 2 i sin ( Q/2 ) + - cos(Q/2) 


^ r ^ [ exp ( i o ) - 1 ] ♦ 


r ^ [ 1 -exp( -i* ) ] 

l 1 C j + 2 Tj 


(9.37) 


- - 2 i sin ( Q/2 ) + j cos(Q/2) ^(r. + 7 ) [exp(i<t>) — 1] -i- 

(r. - 7 77 ) [ 1-exp ( -i$ ) ] I + H.O.T. 

3 2 3 ] ) 


■ -2i sin(Q/2) + ir cos(Q/2) sinU) + 

4 7 ^ cos ( Q /2 ) [ cos ( S ) - 1 ] + H.O.T. 

2 3 j 


(9.38) 


The H.O.T. are neglected because they are of the 
same order of magnitude as other terms already neglected in 
the derivation of the asymptotic amplitude equation. To 
this same level of asymptotic accuracy a 1# a 2 and a, are 
exactly the same as in the analysis of method 1 . 


so now 


motion 


The dispersion relation remains. 


tan ( Q /2 ) » - sin ( $ ) 


(9.26) 


7 q-r cos(Q/2) [cos(« )-1] 


(9.39) 


°° 2 3 j 

Afte:* some algebra it follows that the equations of 


are , 


* r cos(t) cos 1 ( !l /2 ) 
an 


(9.40) 


d® 1 3 r . . . 5 i y \ 

— ■ - — — sinO) cos 2 (Q/2) 
dn 2 3 z 


(9.41) 


4 “ In ( E ) - 7 ^ [ 1 -cos ( 0 ) ] cos 2 ( Q/2 ) 
an d 3 


(9.42) 


Method 3 is a Trapezoidal method which is second 


IOC 
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order accurate in Ax for cases with smoothly varying Ax. 


t « t - 


r . i +r , . 

:+* l-fr 


u A + — — w 7 ] o 

tx r . . , +r . , tx 3 


n + * 


r 3^ +r 3-J- 


(9.43) 


The dispersion relation is the same as for methods 
1 and 2 so after some algebra, 


a, - r-r cos(fl/ 2 ) [cos(»)- 1 ] 
3 3 

and the equations of motion are , 

4 ^ * r cos(«) cos*(Q/ 2 ) 
dn 


d« 3r ... 

— ■ - — sin(®) cos 2 (0/2) 
dn 83 


(9.44) 


4- In ( E ) - 2 — [ 1 *cos ( ♦ ) ] cos 2 ( 0/ 2 ) 
dn 3 j 


(9.45 ) 


(9.46) 


(9.47) 


When the vavepackets of methods 1-3 reach the 
downstream boudary at j »200 they are reflected into 
backward travelling wavepackets. The energy E, and the 

wavenumber ®, of the reflected wavepacket are related to 
the energy E t and wavenumber ® t of the incident wavepacket 
> */ the equations, 

(9.48) 


In ( E 2 ) - In ( Ej ) + 2 In | R | 

U 


where R is the amplitude reflection coefficient. 

u 

For space extrapolation (see 93.5.2) 

sin ( ♦ l / 2 ) 

1 R j' * sin ( ® 2 / 2 ) 


( 9.49) 


tan ( ® j / 2 ) 


(9.50) 


101 


OF POOR QUALITY 


For space-time extrapolation (see S3. 5. 3) 



sin [ ( » x -0 )/2] 
sin [ ( $ 2 -Q ) / 2 ] 


( 9 . 51 ) 


For the box boundary condition (see 53.5.4) 


I R„l 


cos(» x /2) tan(Q/2) - r sin(« 1 /2) 
C03 («j/ 2) tan(Q/2) - r sin(» x /2) 


- tan 3 ( /2 ) (9.52) 

after substituting for tan(fl/2) using the dispersion 


relation and replacing 


The reflection relations at the upstream boundary 

are , 

(9.53) 

In ( E x ) =» ln(E : ) (9.54) 

since R * -1 (see S3. 5.1) 

0 


Method 4 is the Box method discussed in S3. 3. 2 

0 (9.55) 


[ U 5 + r. u a ] 

x t j t x ] + £ 


The dispersion relation is 


tan(Q/2) » r tan(o/2) 


(9.56) 


and the equations of motion for the wavepacket are 
r [ 1 + tan 2 (»/2) ] cos*(Q/2) 


ii 

dn 


~ * -2 tanU/2) cos J (Q/2) 
an o'] 


(9.57) 


(9.58) 


1 n ( E ) » -2 t an 2 ( ♦ / 2 ) cos l (fl/2) 
3 3 
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9.3 Numerical Results 


9.3.1 Trapezoidal Method with Variable CFL Number 


This example uses. 


Method type » 2 ; one of the Trapezoidal methods 

Boundary type - 1 ; space extrapolation 

0.2 r , _ ■ 0.04 


r 0 " °- 05 


200 


cnt 


Figure 6 shows X(n) and ln[E(n)] both predicted and 
experimental. This example shows the movement of a 
wavep.'r.ket and the change in its er.ergy due to the variation 
in the CFL number. The agreement between the predicted and 
experimental values is excellent. The energy of the 
analytic solution id constant so the wavepacket theory has 
successfully predicted almost all of change in the numerical 
energy due to variable Ax. 
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9.3.2 Box Method with Variable CFL Number 


This example uses, 


Method type » 4 j Box method 


r . » 0.05 r_ rtrt * 0.2 
0 200 


crit 


0 .04 


Figure 7 shows X(n) and ln[E(n)]. 
agreement between experiment and theory is 


As in 99.3. 
excel lent . 


/ 


1 the 


. 50 
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9.3.3 Trapezoidal Metnod with Space Extrapolation 


This example uses, 


Method type ■ 1 j one of the Trapezoidal methods 

Boundary type ■ 1 ; space extrapolation 


r 0 “ 1 ’° 


200 


1 . 0 


r . „ - 0 .3 
cri t 


Figure 8 shews X(n) and ln[E(n)] . This example 
illustrates the effect of the downstream boundary reflecting 
a wavepacket with reduced energy. Because of the finite 
length of the wavepacket the drop in energy is smeared and 
X(n) does not quite reach 1.0 . 
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9.3.4 Trapezoidal Method with Space-time Extrapolat i o_n 

This axampls is the seme as 99.3.3 except that 

Boundary type » 2 j space-ti?.« extrapolation 

Figure 9 shows X(n) end lr. [E(ni] . The energy of 
the reflected wavepacket is less than in 99.3.3. As a 
consequence the first order terms which are neglected in 
the asymptotic boundary conditions are more significant and 
so the energy decrease is more smeared and there is a 
larger discrepancy between the experimental and predicted 
values . 


• / i \ 







9.3.5 Trapezoidal Method with Box Boundary Condition 


This example is the same as 59.3.3 except that 


Boundary type * 3 ; box boundary condition 


Figure 10 shows X(n) and in[E(n)]. The energy 
in this example is three times that in 59.3.3 because, 


In ( E j ) - In ( E x ) + 2 ln( | R | ) 

J 


and 


I tan(»,/2) space extrapolation 

I tan 1 (♦ l /2) space-time extrapolation 


Thus the box boundary condition increases the 


drop 

(9.49) 

(9.50) 
(9.52) 


overall convergence rate by factor 3 with minimal 
computational effort. 
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3C0 


0 







9.3.6 Wavepacket Outflow in Box Method 


This example uses, 


Method type 
r 0 * °* 4 


■ 4 

r 2 0 0 


; Box method 

0.4 r * 0.2 

cr it 


Figure 11 shows X(n) and ln[E(n)] . Note that 
the wavepacket reaches the downstream boundary the 
experimental value for X(n) remains near 1.0 and lr.(E) 
decreases rapidly towards . 


when 



•Oil 








9.3.7 Instability of Trapezoidal Method with Space-time 

Extraoolation 


Thi3 example uses, 


Method type 
Boundary, type 

r 0 " 3 *° 


- 1 

- 2 

r 2 0 0 


; one of the Trapezoidal methods 
; space-time extrapolation 

-3.0 r - 2.1 

crit 


Figure 12 show3 X(n) and ln[E(n)]. Both the 
theoretical prediction and the numerical result show that 
the energy increases every time the wavepacket reflects off 
the downstream boundary and so the numerical scheme is 
unstable. The CFL stability condition for the case with 
constant CFL number is obtained by considering the 
amplitude reflection coefficient. 


|R j' 



sin 


l 2 , 


(9.52) 


The dispersion relation is 


tan(3/2) - | sin(«) (9.26) 

Now consider the two cases r<2 and r>2 . 


a ) r < 2 

0 < » : < */2 --> tan(3/2) < 1 

--> 3 < ir / 2 

--> ( t / 2 - 3) > 0 and (*/2 - $ L ) > 0 


--> ir > | ( ir/2-3 ) + 



(ir/2-« l ) | > | ( ir / 2 




> 

sin — 


l 2 J 


3 ) - ( ir/2-o l ) | > 0 


I 



PvjJ 

1 0 ICO 


FIGURE 12. EXAMPLE 9. 
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3.7 TRAPE 


SPACE-TIME 







ONIGINAL PACE rs 
0F POOR QUALITY 


n > 
n < 9 

As 9 2 varies from ir/2 to v , 
continuously so at some intermediate 

£1 - 9 2 

and hence | R | ■ ® 

J 

The numerical scheme is thus stable if, and only if, 
the CFL number is less than 2. This stability condition 
has previously been derived by Beam, Warming and Yee [5] . 


ir/2 ««> £1>9 

Q-o, varies 
value , 
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— > i R ji < 1 

b) r > 2 

For 9 * if / 2 , tan(Q/2) > 1 **> 

For 9“ff , tan(Q/2) ■ 0 *»> 


In a numerical experiment an infinite amplitude 
reflection coefficient does not occur because the first 
order derivatives of the amplitude which are neglected in 
the asymptotic boundary conditions become significant. In 
fact in all the unstable cases I have tried the agreement 
between experiment and prediction is poor because of the 
neglected first order terms in the boundary condition and 
the neglected second order terms in the amplitude equation. 
The example given is ono of the best. The qualitative 
effects of these neglected terms can be understood as 
follows; 


_ . 3 A 

Since — 
3 n 


amplitude equation 
3 2 A 


3 A 

r — all the second order terms in the 
g 3 j 


3 2 A 
3 n 2 


can be expressed in 


terms of 


3j 2 


3 2 A 3 2 A 

3 n3 j ' 3 j 2 

Thus the amplitude equation including 


second order terms has the form, 


v 


3 2 A 

3 j 2 


3 A 3 A 

- — + r — 
3 n g 3 j 


(9.61) 


where v is a function of fl , $ and corresponds to an 
artificial viscosity. The effect of this artificial 
viscosity is to smear the wavepacket increasing its length 
and decreasing its maximum amplitude. This has the largest 
effect on X(n) since the longer the wavepacket the further 
X(n) must be from the ends of the domain. X(n) still 
oscillates approximately in phase with the predictions but 
the amplitude of the oscillations decreases steadily. The 
effect on t^e energy is much smaller. 

The downstream boundary condition including first 
order terms can be written as, 

A 2 ( j , n ) + t, |^ l (J,n) - r [A^J.n) + T t ^ l (J,n)] (9.62) 

d R J dn 

where r^tj are functions of Q , t . If and T a are both 

small compared to T then (9.62) is approximately equal to 

A 

A, ( J,n+r, ) * R A l (j,n+r l ) (9.63) 

U 

Thus the amplitude is reflected with a delay of 
Tj-Tj . This explains the fact that in almost all the 
examples in this chapter the reflected wavepacket lags 
behind the position predicted by wavepacket theory. 
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9.3.8 Instability of Trapezoidal Method Due to Varying Mesh 


This example uses. 


Method type » 3 
Boundary type * 1 


r o ■ ’•° 


r 200 


; one of the Trapezoidal methods 
; space extrapolation 

-10.0 r . -0.8 

cr it 


Figure 13 shows X(n) and ln[E(n)] . The agreement 
between experiment and prediction is good for the energy 
but as in 59.3.7 the agreement is poor for X(n) because of 
the effect of the second order terms which are neglected in 
the asymptotic amplitude equation. The significance of 
this example is that this numerical scheme is stable for 
uniform meshes which give constant C p L number r but if the 
mesh, and hence r, varies sufficiently as in this example 
the scheme becomes unstable. This instability is best 
understood by expanding the finite differ^ uvO equation in 
computational space. 


5 + 

t r . 




+ r 


j+r j 




y A ♦ 
t x 


M ■ 


r j** 


+ r 


.1 




(9.43) 


r . 


so [ «t * T 


3 r . n + i- 

U + — u 6 1 * H.O.T. ] U. 2 - 0 

t 2x 3 ] t x j 


(9.64) 


3 r , 

The term — u 5 Z corresponds to a viscous term in 
3 j t x 

3 r 

computational space. If — is positive it corresponds to 
negative artificial viscosity and so causes the instability 
in the above example. 
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9.4 Program Listing 


The program is written in FORTRAN 4 PLUS to be run on a 
PDP 11-70 with Versatec graphics subroutines. 

C ***** PROGRAM WAVE 

DIMENSION U( 0 : 202 ) ,V( 0:202 ) ,A( 0:200,3 ) ,W(0: 200) ,XJ( 0:20 2) 
DIMENSION EX (402) , EE (402) f TX(402) , TE( 402 ) ,T( 402 ) 

EQUIVALENCE ( TX ( 1 ) ,U( 0 ) ) , ( TX ( 204 ) , V{ 0 ) ) , (TE( 1 ) , A( 0 , 1 ) ) , 

1 (T( 1 ) ,A(0,3) , (T(202) ,W(0) ) 


EXTERNAL R,J 
REAL K, J 

COMMON /RCONST/RO , R1 ,MT 
COMMON /X/X( 0:200) 


C ***** Input parameters 

TYPE * , ' INPUT TERMINAL TYPE' 

TYPE *,'0 VERSATEC' 

TYPE *,'3 VT100 WITH GRAPHICS' 

TYPE * , ' 4 CHROMATICS ' 

TYPE *, ' 5 VISUAL 500 ' 

ACCEPT * , NT 

TYPE *, ' ' 

TYPE *, 'INPUT METHOD TYPE' 

TYPE *,'1-3 DIFFERENT TRAPEZOIDAL METHODS - SEE NODES' 
TYPE * , ’ 4 BOX METHOD ' 

ACCEPT * , MT 

IF (MT.EQ.4) GOTO 1 
TYPE *, ' ' 

TYPE *, 'INPUT DOWNSTREAM BOUNDARY TYPE' 

TYPE *,'1 SPACE EXTRAPOLATION' 

TYPE *,'2 SPACE-TIME EXTRAPOLATION' 

TYPE * , ' 3 BOX METHOD ' 

ACCEPT * , MDB 

1 TYPE * , ' ' 

TYPE *, 'INPUT CFL NUMBERS R( 0 ) , R( 200 ) , RCRIT’ 

ACCEPT -*,R0,R1 , RC 


C ***** Omega definition explained in notes 


OMEGA- 2 . * ATAN ( 0 . 5 * RC ) 


PI-3.14159 
JMAX-200 
M— 400 
XJ( 0)-0. 

DO 2 J1— 1 , JMAX 
DX-1 ./R(FLOAT( J1 ) -0 . 5 ) 
XJ( J1 )-XJ(J1-1 )+DX 
U( J1 )-0. 

V( J1 )-0. 

DO 3 J1- 1 , JMAX 

XJ( J1 )-XJ( J1 )/XJ( JMAX) 

X( JD-XJ(JI) 


Initialise U( j , 0 )+iV( j ,0 ) 

PSI-0. 

J2-JMAX/2 

DO 4 J1-J2-40,J2+40 

IF(MT.LE.3) PHI-ASIN ( RC/R( FLOAT ( J 1 ) ) ) 
IF(MT.EQ.4) PHI-2 . *A7AN ( 0 . 5*RC/R( FLOAT ( J 1 ) ) ) 
PSI-PSI+PHI 

AMP-EXP( -(J1-J2) **2/200. ) 

U(o 1 ) — AMP*COS( PSI ) 

V( J1)-AMP*SIN'PSI) 

W(0)-0. 

KOUNT-O 

TYPE * , ' NO. OF STEPS TILL NEXT PLOT OF U? ' 
ACCEPT * , NSTEP 
IF (NSTEP.LE.0) GOTO 5 

DO 6 KOUNT 2 - 1 , N STEP 
KOUNT- (COUNT-*- 1 
IF (KOJNT.GT.M) GOTO 7 

Calculate new u+iV 

CALL METHOD ( U , A , W , JMAX , R , MT , MDB ) 

CALL METHOD(V,A,W, JMAX, R,MT, MDB) 


Calculate new X,Log(E) of wavepacket 

S-0. 

SJ-0. 

XJ ( JMAX*- 1 ) -XJ ( JMAX ) 
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DO 8 J1-1,JMAX 

P-(U(J1)**2+V( J1)**2)*(XJ(J1+1)-XJ( J1-1) )/2. 

S-S+P 

8 SJ«SJ+P*XJ(J1) 

EX ( KOUNT )-SJ/S 

6 EE( KOUNT )-LOG(S) 

5 CALL OUTPTKXJ ,203, 'X ','U *,NT) 

GOTO 9 

C ***** Calculate predicted X( n ) ,Log( E( n ) ) 

C ***** Initial valu j(1) is passed to prediction subroutine PRED 

C ***** as TXM ) 

7 TX(1)-J(EX( 1) ) 

TE( 1 )-EE( 1 ) 

CALL PRED ( T , TX , TE , JMAX , M , RC , R , MT , MDB ' 

C ***** Plot results 

12 TYPE *, ’PLOT E,X7 ' 

ACCEPT 1000 , C 
IF1C.EQ. 'N ' ) GOTO 10 

IF ( C . EQ . ' E ') CALL OUTPT2 ( j.’ , EE ,TE , 402 , ' N ',’LN E ' , NT ) 

IFtC.EQ.'X ') CALL OUTPT2 l T , EX ,TX , 402 , 'N ','X ' ,NT) 

IFfC.EQ. 'Y ' ) GOTO 11 


GOTO 12 

11 CALL OUTPT3 ( T , EX ,TX , EE ,TE, 402 , 'N ’ , 'X ' , ' LN E ' , NT ) 

10 CALL PLOT(0. ,0. ,999) 

STOP 

1000 FORMAT(A4) 

END 


C ***** Function calculates R(j) 

FUNCTION R(J) 

COMMON /RCONST/RO , R1 ,MT 
REAL J 

t 

R-R0*EXP(LOG(R1/R0)*J/200. ) 
RETURN 


* 

•( 


END 
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C ***** Function calculates X ( j ) 

FUNCTION X(J) 

COMMON /X/XJ( 0:200) 

REAL J 

IF( J.GT . 200 . ) J-200. 

IF(J.LT.0.) J-0. 

J1-INT( J) 

X»XJ ( J 1 )+( J- FLOAT ( J1 ) ) * ( XJ( J1+1 ) -XJ( J 1 ) ) 

RETURN 

END 


C ***** Function calculates j(X) 

FUNCTION J(X) 

COMMON /X/XJ( 0:200; 

REAL J 

IF(X.GT.1.) X-1. 

IFtX.LT.O.) X-0. 

J1»0 

1 J1-J1+1 

IF ( X J ( J 1 ) • LT . X ) GOTO 1 

J-FLOAT( J1)-(XJ( J1 )-X)/(XJ( J1 )-XJ(J1-1) ) 

RETURN 

END 


SUBROUTINE METHOD (U,A,W, JMAX , R , MT , MDB ) 

C ***** METHOD sets up the coefficients of the tridiagonal equations 

for the calculation of the new U using method MT and downstream 
boundary type MDB, if needed. The tridiagonal equations are 
C solved by TRID and the new values of U are returned to T. 

DIMENSION U ( 0 : JMAX ) ,A(0:JMAX,3) ,W(0:JMAX) 

IF(MT.EQ.4) GOTO 1 

C ***** Set up coefficients for Trapezoidal and Backward Euler interior 
C ***** schemes 


W I 
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DO 2 J-1 , JMAX-1 
GOTO( 3,4,5) MT 

3 Cl— 0.25*R(FLOAT( J) ) 

C2— Cl 

GOTO 6 

4 Cl— 0.25*R(FLOAT( J)-0.5) 

C2-0.25*R(FLOAT( J) -*-0 . 5 ) 

GOTO 6 

5 R1-R(FLOAT( JJ-0.5) 

R2=»R( FLOAT ( J)*0.5) 

Cl— 0.5*R1**2/(R1+R2) 

C2-0.5*R2**2/(R1+R2) 

6 A(J,1)*C1 

A( J,2)-1 . -C1-C2 
A( J,3)-C2 

2 W(J)»J(J)+C1*(U(J)-U(J-1) )+C2*(U(J)-U(J+1) ) 

GOTO (7,8,1) MDB 


C ***** Set up coefficients for space and space-time extrapolation at 
C downstream boundary 

7 A ( JMAX , 1 ) — 1 . 

A ( JMAX , 2 ) * 1 . 

W(JMAX)«0. 

GOTO 9 

a A( JMAX , 1 ) *0 . 

A ( JMAX , 2 ) “ 1 . 

W( JMAX) «0( JMAX-1 ) 

GOTO 9 


C ***** set up coefficients for box method on interior or at downstream 
C boundary as appropriate 

1 JMIN-JMAX 

IF ( MT . EQ . 4 ) JMIN* 1 
DO 10 J* JMIN , JMAX 
C*R( FLOAT ( J) -0 .5 ) 

A( J, 1 )=»1 .-C 
A(J , 2 ) * 1 . +C 
A(J,3)»0. 

10 W(J)»D(J-1)*A(J,2)+U(J)*A(J,1) 


*** Se up coefficients for upstream boundary 

9 A(0,2)-1. 

A(0,3)»0 . 

W(0)-0. 


C 
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14 CALL TRID ( U , A , W , JMAX ) 

RETURN 

END 


SUBROUTINE TRID ( U , A , W , JMAX ) 

DIMENSION U{ 0 : JMAX) ,A(0:JMAX,3) ,W(0: JMAX) 

DO 1 J— 1 , JMAX 
OA(J,1)/A(J-1,2) 

A(J,2)«A( J,2)-A( J-1 ,3)*C 

1 W( J)»W{ J)-W( J-1 )*C 
U(JMAX)»W( JMAX)/A( JMAX, 2) 

DO 2 J- JMAX- 1 ,0,-1 

2 U(J)«(W(J)-U(J+1)*A(J,3) )/A( J,2) 

RETURN 

END 


SUBROUTINE PRED ( T , TX , TE , JMAX , M , RC , R , MT , MDB ) 

DIMENSION T ( M ) , TX ( M ) , TE ( M ) 

REAL J,K 
EXTERNAL R,X 

DR(J)*100.*(R(J+O.O n 5)-R(J-0.005)) 

PI-3. 14159 
J-TX ( 1 ) 

TX ( 1 )-v( J) 

T( 1 )-1 . 

IF(MT.EQ.4) GOTO 1 


C ***** Prediction for trapezoidal schemes 

X=»ASIN(RC/R(J) ) 

IF(RC.LT.O.) K-PI-K 
OM-ATAN ( RC/2 . ) *2 . 

C1-1./( 1 .+0.25*RC**2) 

DO 2 KOUNT-2 ,M 
T ( KOUNT ) =* FLOAT ( KOUNT ) 

DJ-R( J)*COS(K)*C1 


DK— DR( J)*SIN(K)*C1 ' 

J- J+0 . 5 * ( DJ+ R( J+DJ ) *COS ( K+DK ) +C 1 ) 

K-K+O . 5 * ( DK- DR ( J+DJ ) + blN ( K+PK ) +C 1 ) 

TX ( KOUNT) «X( J) 

DE«DR( J+0.5*DJ)*C1*( 1 . -COS( K-0 .5+DK) ) 

IF(MT.EQ.I) TE ( KOUNT ) — TE ( KOUNT- 1 ) 

IF( MT.EQ . 2 ) TE ( KOUNT )«TE( KOUNT- 1)+DE 
IF(MT.EQ.3) TE ( KOUNT )»TE( KOUNT- 1) +2. *DE 

IF(J.GT.O.) GOTO 3 
J— J 

K-ASIN( RC/R( J) ) 

TX ( KOUNT ) ( J ) 

IF ( J . LT . FLOAT ( JMAX ) ) GOTO 2 
J- 2 . *FLOAT ( JMAX ) - J 
K— PI-ASIN( RC/R( J) ) 

TX ( KOUNT ! -X ( J ) 

GOTO (6,7,8) MDB 

TE( KOUNT )-TE( KOUNT) -2. *LOG( TAN ( K/2 . ) ) 

GOTO 2 

TE ( KOUNT )-TE (KOUNT) -2. *LOG(ABS( SIN ( (K-OM)/2. )/SIN( ( PI-K-OM)/2. ) ) ) 
GOTO 2 

TE ( KOUNT ) =TE ( KOUNT ) -6. *LOG( TAN (K/2. ) ) 


CONTINUE 

RETURN 


Prediction for box scheme 

T1-0.25*RC*RC 
T2*1 ./( 1 .+T1 ) 

T3-0 .5*RC*T2 

K=2 . *A1 AN ( 0 . 5*RC/R( J ) ) 

DO 9 KOUNT- 2, M 
T ( KOUNT I =■ FLOAT ( KOUNT ) 

IF( J.GE.FLOAT( JMAX) ) GOTO 10 
R1=R( J) 

DJ=(R1+T1/R1)*T2 
DK--2 . *DR( J) +T3/R1 
R2-R( J+DJ) 

J-J+0 .5* ( DJ+( R2+T1/R2 ) *T2 ) 

K-K+0 . 5 * ( DK+DR( J+DJ ) *T3/R2 ) 

TX ( KOUNT) «X( J) 

TE( KOUNT) «TE( KOUNT- 1 ) -2. *T1*T2+DR( J+0. 5*DJ ) /{ R1 *R2 ) 
GOTO 9 


V 
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10 TX(KOUNT)*1. 

TE ( KOUNT ) «TE ( KOUNT- 1 ) 

9 CONTINUE 

RETURN 

END 


SUBROUTINE OUTPT1 ( X,Y,NPLUS2 , Cl ,C2,NT) 

DIMENSION X ( NPLUS2 ) , Y ( NPLUS2 ) ,XD(40) 

N-NPLUS2-2 
NPLUS1-N+1 
XL* 5 . 

YL*4 . 

CALL PLOTS ( 0 ,0 ,NT) 

CALL SCALE ( X , XL , N , 1 ) 

CALL SCALE ( Y , YL , N , 1 ) 

DO 1 I* 1 ,N/8 

1 XD( I )-(X( 8*1+1 )-X( 8*1-7) ) /X( NPLUS2 ) 

IF(NT.EQ.O) GOTO 2 
CALL FACTOR (1.8) 

CALL PLOT( 1 . , 1 . , -3 ) 

CALL AXIS( 0 . ,0. ,C1 ,-4,XL,0. ,X(NPLUS1) ,X(NPLUS2) ) 
CALL AXIS( 0. ,0. ,C2,4,YL,90. ,Y(NPLUS1 ) ,Y(NPLUS2) ) 
CALL LINE( X , Y ,N , 1,1,3) 

CALL GRID ( 0. ,0. , 1 000+N/8 ,XD , - 1 , YL , - 1 ) 

GOTO 3 

2 CALL 1ACTORM.4) 

CALL PLOT ( 5 . , 1 . ,-3) 

CALL AXIS( 0 \ ,0. ,C1 , -4, XL, 90 . ,X(NPLUS1 ) ,X(NPLUS2) ) 
CALL AXIS( 0. ,0. ,C2,4,YL, 180. ,Y(NPLUS1 ) ,Y(NPLUS2) ) 
Y ( NPLUS2 ) ■- Y ( NPLUS2 ) 

CALI. NEWPEN ( 2 ) 

CALL LINE (Y,X,N,1,1,3) 

CALL NEWPEN ( 1 ) 

CALL GRID ( 0 . , 0 . , - 1 , - YL , 1 0OO+N/8 ,XD , - 1 ) 

3 CALL PLOT(0. ,0. ,-999) 

RETURN 

END 
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SUBROUTINE OUTPT2 ( X , Y 1 , Y2 , NPLUS2 , C 1 , C2 , NT ) 

DIMENSION X ( NPLUS2 ) , Y 1 ( NPLUS2 ) , Y2 ( NPLUS2 ) , Y ( 4 ) 

N-NPLUS2-2 

NPLUS 1 *N+ 1 

XL=5. 

YL=4. 

Y( D-1.E10 
Y(2) — 1 .E10 
DO 1 J*1,N 

Y( 1 )«AMIN1 (Y( 1 ) ,Y1 (J) ,Y2( J) ) 

1 Y( 2 ) = AMAX 1 { Y ( 2 ) , Y1 ( J) , Y2 ( J) ) 

CALL PLOTS ( 0 ,0 ,NT) 

CALL FACTO R( 2.0) 

CALL PLOT( 1 . ,0.75, -3) 

CALL NEWPEN ( 2 ) 

CALL SCALE ( X , XL , N , 1 ) 

CALL SCALE ( Y , YL ,2,1) 

CALL AX I S ( 0 . ,0. ,C1 ,-4, XL , 0 . ,X( NPLUS 1 ) ,X(NPLUS2) ) 
CALL AXIS ( 0 . , 0 . , C2 , 4 , YL ,90. ,Y(3) ,Y(4) ) 

Y 1 ( NPLUS 1 ) =Y ( 3 ) 

Y1 (NPLUS2)=*Y(4) 

Y2 ( NPLUS 1 ) “Y( 3 ) 

Y2 ( NPLUS2 ) =* Y ( 4 ) 

CALL NEWPEN ( 3 ) 

CALL LINE(X,Y1,N, 1,25,1) 

CALL LINE(X,Y2,N, 1,0,0) 

CALL PLOT(0. ,0. ,-999) 

RETURN 

END 


SUBROUTINE OUTPT3(X,Y1 ,Y2,Z1 , Z2 ,NPLUS2 , CX ,CY ,CZ ,NT) 

DIMENSION X(NPLUS2) ,Y1(NPLUS2) ,Y2(NPLUS2) ,Y(4) 
DIMENSION ZKNPLUS2) , Z2 ( NPLUS2 ) , Z ( 4 ) 

EQUIVALENCE (Y(1),Z(1)) 

N-NPLUS2-2 
NPLUS 1=N+1 
XL»4. 

YL=2. 

ZL=2 . 


( 


• ■ I 

OF POOR QUALITY 


CALL PLOTS ( 0 , 0 , NT ) 

CALL FACTOR (1.9) 

CALL SCALE ( X , XL , N , 1 ) 

Y( 1 )— 1 -E10 
Y(2) — 1.E10 
DO 1 J-1,N 

Y( 1 J-AMIN1 (Y( 1) ,Y1 (J) ,Y2(J) ) 

Y(2)»AMAX1 (Y(2) ,Y1 ( J) ,Y2 ( J) ) 

CALL PLOT (3. 5, 0.5, -3) 

CALL SCALE ( Y , YL ,2,1) 

CALL NEWPEN ( 2 ) 

CALL AXIS(0. ,0. ,CX , -4 ,XL, 90 . ,X( NPLUS 1 ) ,X(NPLUS2) ) 

CALL AXIS( 0 . , 0 . ,CY,4,YL, 180. ,Y( 3) ,Y(4) ) 

Y1 (NPLUS1 )-Y(3) 

YMNPLUS2)— Y(4) 

Y2 ( NPLUS 1 ) *Y ( 3 ) 

Y2 (NPLUS2 ) «-Y( 4 ) 

CALL MEWPEN(3) 

CALL LINE( Y1,X,N, 1,25,1) 

CALL LINE (Y2,X,N, 1,0,0) 

CALL ’EWPEN(I) 

CALL v J r T O(0. ,0. ,2,-1. ,4, 1 . ,-21846) 

CALL NEWPEN ( 2 ) 

CALL PLOT (-2. 3, 0.4, 3) 

CALL PLOT(-2.3, 1.0,2) 

CALL SYMBOL(-2.23, 1 . 1 ,0. 14, 'WAVEPACKET THEORY 90 ., 17 ) 
CALL PLOT (-2. 6, 0.4, 3) 

CALL SYMBOL ( -2. 6, 0.5, 0.08, 1,0., -2) 

CALL SYMBOL (-2.6,0.7,0.08,1,0. ,-2) 

CALL SYMBOL (-2.6,0.9,0.08,1 ,0. ,-2) 

CALL PLOT( -2.6, 1 .0,2) 

CALL SYMBOL(-2.53, 1 . 1 ,0. 14, 'NUMERICAL EXPERIMENT 90 ., 20 ) 

Z( 1 )«1 .E10 
Z(2) — 1.E10 
DO 2 J-1 ,N 

Z( 1 )-AMIN1 (Z( 1 ) , Z 1 ( J ) , Z2 ( J ) ) 

Z( 2)»AMAX1 (Z(2) ,Z1( J) ,Z2( J) ) 

CALL PLOT (3.2,0. , -3 ) 

CALL SCALE ( Z , ZL , 2 , 1 ) 

CALL AXIS(0. ,0. ,CX, -4 ,XL , 90 . ,X(NPLUS1 ) ,X(NPLUS2) ) 

CALL AXIS( 0. ,0. , CZ , 4 , ZL , 180. ,Z(3) ,Z(4) ) 

ZKNPLUS1 )-Z(3) 

Z1 (NPLUS2 ) *-Z( 4 ) 

Z2 ( NPLUS 1 ) *Z ( 3 ) 

Z2 ( NPLUS2 ) »- Z ( 4 ) 

CALL NEWPEN ( 3 ) 

CALL LINE(Z1 ,X,N, 1 ,25, 1 ) 
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CALL LINE(Z2,X,N, 1,0,0) 

CALL NEWPEN ( 1 ) 

CALL GRID ( 0 . ,0. ,2,-1. ,4,1. ,-21846) 

RETURN 

END 
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10. Conclusions 


The validity of the asymptotic approach developed in 
this paper is demonstrated by the numerical results in 
chapter 9. The limitations of the wavepacket theory are due 
to the asymptotic approximations involved in treating the 
wavepacket as a particle. The stability analyses in 
chapters 5 and 6 use fewer approximations and so the 
asymptotic errors will be substantially smaller. In 
particular when the coefficients are constant the analysis 
in chapter 6 reduces to the P-stability analysis of Beam, 
Warming and Yee [5] . 

The calculation of the asymptotic amplitude equation 
and asymptotic boundary conditions for a particular case is 
no more difficult than a normal Von Neumann analysis. For 
applicable cases the wavepacket theory and the stability 
analysis of chapter 5 are straightforward. The general 
stability analysis of chapter 6 will usually require 
numerical computation. In the more complex cases the main 
benefit from this theory will be the insight given by the 
asymptotic amplitude equation and boundary conditions. The 
amplitude equation gives the group velocities of the 
different wa/enumbers and the effect of varying 
coefficients, which is of great interest since in 2-D 
cascade geometries cell lengths can vary by factors of up to 
100 in inviscid calculations and 1000 in viscous 
calculations. The asymptotic boundary conditions give the 
amplitude reflection coefficients which provide a practical 
criterion for choosing the best numerical boundary 
conditions . 


There are various possibilities for future research 


in this area. Further applications to relatively simple 
problems can be done to gain insight into understanding 
harder problems and improving boundary conditions. 

Numerical procedures, such as those suggested in 
§6. 1.1-6. 1.3 , can be developed to solve the equations given 

by the stability analysis in chapter 6. Finally tre 
asymptotic amplitude equation and boundary conditions can be 
extended to 2-D and 3-D. 
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Appendix 


A . 1 Finite Difference Operator Notation 


An operator notation for finite difference equations 
simplifies the analysis of finite difference schemes and is 
a necessity for making any general statements and proving 
them . 

The principal operators are central difference, 

central averaging E^ , shift operator, A ^ , forward 

difference and 7 , backward difference. Their definitions 

x 

are; 


5 0, ■ a, - U. 

mx j j+ra/2 j-m/2 


( A . 1 ( a ) ) 


U mx U j " 2 (U j+m/2 * U j-m/2 ) 


E 0, ■ U, 
mx j j+m 


a a . - a . - u . 

ax ] j+m j 


7 U a - - 0. 

mx j j j-m 


( A. 1 ( b) ) 
(A. 1 ( C) ) 
(A. 1 ( d) ) 
( A. 1 ( e ) ) 


Usually these definitions will be used with m» 1 . 

The main exceptions are 6^ x which is a node -ce n t e r e d central 

difference, 


6 _ U a - 0, , - U. , 
2 x j 3 + 1 j-1 


(A. 2) 


and E which can be used to define a general linear 
m x 


operator , 



ORIGINAL PAGE 19 
OF POOR QUALITY 


so 


M 

L . - ) a . E 

j </_, jra mx 

m* 1 


(A. 4) 


When there are several independent variables the 
subscript on the finite operator denotes the direction of 
the shift, differencing or averaging. For example if, 


U . - u( x , t ) 
J J n 


(A. 5 ) 


then 6 0* - 0 n 

x j + f j+1 


- u 


(A. 6) 


and 5 t o"** - o"* 1 - u” 


(A. 7) 


The general shift operator expression for a finite 
operator in 2-D is 


L 


j 


) C ( j ) E E 

rap mx pt 

m , p 


(A. 8) 


In applications however this expression can be very 
complicated and it is usually simpler to express L as a 
polynomial in the finite operators. As an example the 
operator in 53.3.2 has the polynomial form 


l. - u a_ ♦ r, u a 
j x t j t x 


but in the shift operator form it is, 


(A. 9) 


1 + r 

E E 

2 x/2 t/2 


♦ 


1 - r 
2 


E -x/2 E t/2 



E . E 
x/2 -t/2 


lt_£ 

2 


E E 

-x/2 -t/2 


( A . 1 0 ) 


Part of the advantage in using operator notation 
when analysing finite difference schemes arises because all 
of the finite operators have the same eigenfunction which in 
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OmGWAL PAGE 13 
OF POOR QUALIFY 


4 exp (i ( j d-nfl ) ] 

X 

u x exp [ i ( j 9 -na ) ] 
E x exp t i ( j »-nfl ) ] 
exp [ i ( j ♦-nQ ) ] 
7 exp [ i ( j »-nfl ) ] 

and similarly, 

4 1 exp [ i v * -nQ ) ] 

u t exp [ i ( j *-nQ ) ] 

exp r i ( j * -nQ ) ] 

A t exp [ i ( j ♦-nQ ) ] 

7 t exp [ i ( j »-nQ ) ] 


exp [ i ( ( j+fc ) e-nQ ) ] - exp [ i ( ( j -*• ) *-nQ ) ] 


21 sin(t/2) exp [ i ( j ♦ -nQ ) ] (A. 11(a)) 

t 

exp { i ( ( j*}- ) ♦-nQ ) ] + exp ( i ( ( j -$• ) $-nQ ) ] 

(A. 1 1 (b) ) 


coe(t/2) exp [ 1 ( j 9 -nQ ) ] 
exp [ i ( ( j ♦ 1 ) 9 -nQ ) ] 


exp(i«) exp [ i { j 9 - nQ ) ] (A. 11(c)) 

exp [ i ( ( j ♦ 1 ) « -nQ ) ] - exp [i( j 9 -nfl ) ] 

{ exp ( i * ) - 1 > exp [ i ( j 9 -nQ ) ] (A. 11(d)) 

exp [i ( je-nQ ) ] - exp [ 1 ( ( j- 1 ) 9 -nQ ) ] 

{ 1 - exp ( - i 9 ) ) exp Il( j#-nQ) ] (A. 11(e)) 


-2i sin(Q/2) exp [ i ( j * -nQ 5 ] 
cos(Q/2) exp { i ( j • -nQ ) ] 
exp(-iQ) exp ( i ( j 9 -nQ ) ] 

{ exp(-iQ) - 1 ) exp ( i ( j 9 -nQ ) ] 
{ 1 - exp(iC) ) exp [ i ( j 9 -nQ ) ] 


( A. 1 1 ( f ) ) 
( A. 1 1 ( g) ) 
( A. 1 1 (h) ) 
( A. 1 1 ( i ) ) 
( A. 1 1 ( j ) ) 


A. 2 Principle of the Argument 


Let f(z) be an analytic complex function with 
simple poles in a region of the complex z-plane and let C 
be a closed curve in the region. Then the number of zeros 
of F minus the number of poles of f lying inside C is equal 
to , 


1 f f ' ( z) 

I 0 z 

2 *i J f ( z ) 

C 

’ ITT 1 ln ' £ > >e 

- [ ar,(f) i c 

[ ] £ denotes the change as z goes round C 

anticlockwise . 


( A. 12 ) 


arg(f) is defined by, 

f(z) * R exp(i9) R , 9 real R> 0 (A. 13) 

arg( f ) = 9 ( A. 14 ) 

with the restriction that 9 must vary continuously as z goes 
round C. 

The {.roof is given in many standard texts on complex 
analysis, e.g. [6] . This provides a very simple test when 
considering stability problems in which it is sufficient to 
know whether there ara any zeros in a critical region 
without knowing their exart position. This is the basis of 
the Myquist criterion in control theory stability analysis. 
The test can also be performed numerically relatively 
easily. The step size Az in going round C is decreased, if 
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necessary, until |Asrg(f)| < e. Since ~ [ arg(f) ] is 

an integer there is no rounding error. The only possible 
error is if the magnitude of Aarg(f) over one step lies in 
the range 2nir-e < Aarg(f) < 2nw + e for some integer n 

other than zero. Decreasing t reduces the chance of an 
error at the expense of increased computation. e-ir/6 should 
be adequate in most cases. 
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A. 3 Definitions of Norma and Stability 


The np;as used ir* this paper are generalised L, 
norms. For a continuous function u(x,t) defined on 0<x<X 
the norm | |u(t) | | is defined by, 

X 

| | u ( t ) | | 2 - / |u(x,t)|’ a ( x ) dx (A. 15) 


where a ( x ) is a positive non-zero function. 


For a discrete function U. defined on 0<j<J the 

3 

norm | | U n | | which is a function of n is defined by, 


o n l | * - V |o"|> 


a . 


(A. 16 ) 


where a. is a positive non-zero function. 

3 

The stability used in this paper is asymptotic 
Liapounov stability which is defined as. 


Given <5 > 0 there exists e>0 such that 

I I u( 0 ) | | < e =■»> a) | | u< t) | | < 6 

and b) ||u(t)|| ♦ 0 as t-® 


Condition a) is Liapounov stability which limits 
how large an initially small disturbance can become. 
Condition b) is asymptotic stability which specifies that a 
sufficiently small initial disturbance must ultimately tend 
to zero. 


For linear systems of equations an equivalent 
definition is 


a) There exists M>0 such that ||u(t)|| < M ||u(0)|| 
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and b) ||u(t)|| *0 as t-*® 

The corresponding definitions for a discrete 
function ara. 

Given 4>0 there exists e>0 such that 

I |D°I | < c — > «' | |0 B | | < J 

and h ' | |u n | 1 ♦ 0 as n-»® 

and for linear discrete systems, 

a ^ There exists M>9 such that ||u n j| < M ||U°|| 
and b) | | U n | | ♦ 0 as 


